Three-dimensional numerical modeling of fullspace transient electromagnetic responses of water in goaf*

Chang Jiang-Hao, Yu Jing-Cun♦, and Liu Zhi-Xin

Three-dimensional numerical modeling of fullspace transient electromagnetic responses of water in goaf*

Chang Jiang-Hao1, Yu Jing-Cun♦1, and Liu Zhi-Xin1

The full-space transient electromagnetic response of water-filled goaves in coal mines were numerically modeled. Traditional numerical modeling methods cannot be used to simulate the underground full-space transient electromagnetic field. We used multiple transmitting loops instead of the traditional single transmitting loop to load the transmitting loop into Cartesian grids. We improved the method for calculating the z-component of the magnetic field based on the characteristics of full space. Then, we established the fullspace 3D geoelectrical model using geological data for coalmines. In addition, the transient electromagnetic responses of water-filled goaves of variable shape at different locations were simulated by using the finite-difference time-domain (FDTD) method. Moreover, we evaluated the apparent resistivity results. The numerical modeling results suggested that the resistivity differences between the coal seam and its roof and floor greatly affect the distribution of apparent resistivity, resulting in nearly circular contours with the roadway head at the center. The actual distribution of apparent resistivity for different geoelectrical models of water in goaves was consistent with the models. However, when the goaf water was located in one side, a false low-resistivity anomaly would appear on the other side owing to the full-space effect but the response was much weaker. Finally, the modeling results were subsequently confirmed by drilling, suggesting that the proposed method was effective.

Goaf, water, mine transient electromagnetic method, fullspace, finite-difference time-domain method

Introduction

Water in goaves is responsible for many coal mine disasters as it is very difficult to control. Once dug, highpressure water-filled goaves are a safety threat in coal mines. Presently, the transient electromagnetic method is widely used in detecting water-bearing structures in coalmines (Yu, 2007). Nevertheless, the transient electromagnetic response is affected by the full space;hence, the full-space transient electromagnetic field is different from the half-space transient electromagnetic field and, thus, it is necessary to study the full-space transient electromagnetic responses in water-filled goaves. The 3D forward modeling of the transient electromagnetic field is mainly based on the finitedifference method (Wang and Hohmann, 1993), the integral equation method (SanFilipo and Hohmann, 1985; Zhdanov et al., 2006), and the finite-element method (Goldman et al., 1986). The finite-difference time-domain (FDTD) method directly solves the Maxwell equations in the time domain and avoids the complexity of first solving the responses in the frequency domain and then in the time domain. This is done by using the inverse Fourier transform (Sun et al., 2013). Using the 2D Du Fort–Frankel finite-difference method, Oristaglio and Hohmann (1984) and Yan et al. (2002) discussed the diffusion of the transient electromagnetic fields with time. Liu et al. (2013) performed numerical modeling of the transient electromagnetic responses of 2D line sources using the Du Fort–Frankel finitedifference methods for anomalous fields. Using a modified version of the Du Fort–Frankel 3D finitedifference method, Wang and Hohmann (1993) directly solved the half-space transient electromagnetic responses in the time domain and performed 3D forward modeling by using the FDTD method. All the above-mentioned studies adopted analytical solutions for transmitting sources in uniform half-space, whereas Sun et al. (2013) added the current density to the magnetic curl equations. For the roadways in coalmines, Jiang (2008) studied the transient electromagnetic responses of faults and collapsed columns using the 2.5D finite-difference method. Yang and Yue (2008) studied the effect of roadways on the full-space transient electromagnetic method (TEM) using the 3D difference in the magneticfield diffusion equation. The continuous conditions of the electromagnetic field in complex abrupt interfaces are not considered in existing 3D forward modeling methods. The basic approach of the full-space transient electromagnetic numerical modeling is the same as that of half-space, whereas the transmitting loop for fullspace roadway detection, the multidirectional detection, and the full-space boundary conditions are different from the half-space.

Fig.1 The MTEM for advanced detection.

To solve these problems, we propose loading the minor multidirectional transmitting loop into a Cartesian grid and modifying the finite-difference equation according to the full-space characteristics. We build a 3D full-space geoelectrical model by using the FDTD method and study the full-space transient electromagnetic responses of a 3D geoelectrical body, and analyzed the transient electromagnetic responses of water in goaves.

Advanced MTEM

The mine transient electromagnetic method (MTEM) is a time-domain electromagnetic method based on the electromagnetic induction principle. It uses an ungrounded loop to create a primary pulse magnetic field, which then induces a current and a timedependent electromagnetic field. At the pulse intervals, the induced secondary field is received by the loop in fullspace. Subsequently, data processing and analysis of the received signal allows for the characterization of the geological conditions around the roadway. Unlike ground TEM, both the transmitting and receiving in the MTEM are carried out around the roadways of coalmines. A multiturn coincident loop with length less than 3 m is commonly used in coal mines, as shown inFigure1 (Yu, 2007). A coincident loop with length less than 3 m and different number of turns of transmitter and receiver loops is commonly used for multidirectional detection in coal mines.

The common approach to detect water in goaves is to place the transmitter and receiver loops at the heading face of roadways. Owing to the space limitations of the roadways in coalmines, the commonly used array detection technology of tunnel TEM cannot be used in the MTEM. Thus, only the sector detection method is used. The common arrangement of this method is shown in Figure 1, where it can be seen that all the measurement points are located at the heading face. The arrows indicate the normal direction of the transmitting and receiving antennas at each measurement point. In detection, information about the advanced space can be obtained by rotating the transmitter and receiver loops to collect data at multiple angles.

The FDTD method for electromagnetic fields

Time-domain finite-difference equations and difference scheme

Under quasi-static conditions and for a source-free, isotropic, and nonmagnetic medium, the Maxwell’s equations (Kaufman and Keller, 1983) are

where B is the magnetic induction, H is the magnetic field intensity, E is the electric field intensity, σ is the conductivity of the medium, and J is the conduction current density.

Because the displacement current is ignored, equation (2) does not contain the derivative of the electric field with respect to time and, as a result, the explicit difference scheme cannot be established. Thus, Wang and Hohmann (1993) modified equation (2) to

where the first term on the left-hand side of the equation represents the fictitious displacement current. Introducing this term and assigning proper values to γ can lower the requirements for the time steps in the finite-difference iterations.

Equation (3) suggests that only two of the three components of the magnetic induction are independent. The third can be obtained from the two independent components. In the iterations, the x- and y-components of the magnetic induction can be calculated by using equation (1) and the z-component of the magnetic induction can be obtained by substituting the above two components into equation (3) (Wang and Hohmann, 1993).

In the 3D FDTD method, the earth model needs to be discretized. Unlike ground modeling, a full-space model is needed in underground modeling. The fullspace earth model can be divided into nonuniform grids, as shown in Figure 2, where the source is at the center of the model and the grid spacing increases laterally and vertically away from the source. The Yee grid, shown in Figure 3, is used in each discretized grid. The electric field components are placed at the centers of the prism edges with their directions parallel to the prism edges,and the magnetic field components are placed at the centers of the prism faces with their directions parallel to the normal directions of the faces (Yee, 1966). In Figure 3, each magnetic vector is surrounded by four electric fields, forming the curl of the magnetic field. By simulating Faraday’s law of electromagnetic induction, each electric field component is surrounded by four magnetic fields.

Fig.2 Discretization of a 3D full-space earth model (Yee, 1966).

Fig.3 Position of each field component in the Yee grid (Yee, 1966).

The equations for the electric and magnetic fields, and the specific discrete process and difference equation can be obtained by using the component forms in equations (1) and (5) and the central difference to perform an approximate discretization (Wang and Hohmann, 1993).

With respect to the discretization of Bz, the discrete approach, which is different from that in half-space, needs to be used based on the full-space models to maintain the symmetry of the fields in full space. First, the following equation is obtained from equation (1)

Then, from equation (3), we obtain

In the full-space model in Figure 2, the source is located at the center of the model, which is often taken as the origin of the coordinates. The Bzof each unit on the plane where the origin is located—namely, the plane z = 0—is first calculated by using equation (6). Then, the Bzat other positions is calculated by using equation (7). By using the central difference, equation (6) is discretized as follows:

Equation (7) is discretized as

In the iterations, the Bzof each unit on the plane where the source is located (plane z = 0) is first calculated by using equation (8). Then, the Bzof the units located at z ＜ 0 and z ＞ 0 are calculated using equations (9) and (10). This discrete approach ensures that the fields are symmetric about the plane z = 0. In addition, not only it adds the constraint of equation (3) to the iteration process but also makes Bzindependent of the boundary conditions. Therefore, it improves the precision.

Stability conditions and numerical dispersion

The stability conditions must be satisfied in the finitedifference iterations. The stability condition for 3D iterationsis (Yee, 1966; Wang and Hohmann, 1993)

where Δminis the minimum grid spacing. Equation (11) shows that γ is related to the time step and minimum grid spacing. To ensure the field’s diffusivity, the time step shall satisfy the following (Oristaglio and Hohmann, 1984)

In the calculations, the appropriate time space can be obtained first by using equation (12) and then γ to satisfy the stability condition.

Numerical dispersion occurs in the finite-difference numerical modeling, and the condition commonly used to suppress the numerical dispersion is (Alford et al., 1974)

where λminis the minimum wavelength and N is the smallest number of grid points per wavelength. For second-order differences, N shall be not less than 10. Transient electromagnetic fields are wideband fields. However, high-frequency electromagnetic waves are rapidly absorbed by lossy dielectric media; thus, the remaining electromagnetic waves are mainly lowfrequency waves with frequency generally less than 1 MHz, i.e., the wavelength is greater than 300 m. Therefore, the grid size shall be not more than 30 m, and smaller spatial grids are needed for the shallow parts. The transient electromagnetic field is highly variable initially and numerical dispersion occurs during this period. The numerical dispersion error can be reduced by increasing the spatial sampling rate. Figure 4 shows the relative errors for grids of different sizes versus the number of 500 iterations after the current turn-off. It can be seen that the denser the grids, the smaller the error. When the size of grids is 1 m, the error after 500 iterations is less than 0.1%, which satisfies the precision requirements.

Fig.4 Relative error vs. number of iterations for grids of different size.

Initial source and boundary conditions

The common approach for setting up the transmitting source is to use the electromagnetic field generated by a magnetic dipole or electrical dipole source as the initial source. In practice, this approach requires the medium around the transmitting source to be uniform dielectric and the initial time shall be neither too large nor too small. In principle, the initial time shall be sufficiently small for the solution to be effective in the homogeneous full space. However, a too small initial time leads to insufficient spatial sampling rate, and numerical dispersion occurs easily in the shallow parts. Generally, the resistivity distribution around the source is nonuniform because the transmitting loop used by MTEM is on the roadways of coalmines, and thus there is some difference between the transient field distributions generated by the magnetic dipole source and the minor loop source, respectively. Replacing the minor loop source with the magnetic dipole source produces errors; therefore, we used a minor loop source, with the same size as the actual transmitting loop, as the excitation source. To ensure the field’s stability early at the finite-difference iterations, Sun et al. (2013) proposed adding a trapezoidal transmitting current to the loop. This solution is not model-dependent and can be used to simulate the transient electromagnetic responses of any transmitting current.

In the source region, equation (5) can be modified to (Sun et al., 2013)

where Jsis the density of the source current. When applying the loop source, the spatial position of the loop source can be forced to coincide with that of the electricfield component.

The expressions for the component of equation (14) are

Equation (17) is different from the corresponding equation obtained by Sun et al. (2013). This is because the transmitting loop used in ground TEM is generally placed on the ground and its transmitting current has only horizontal components, whereas the transmitting loop used in MTEM needs to be placed in different directions according to the actual conditions. In the source region, the transmitting current and the iterationsfor Ezneed to be considered.

The difference equations of Exand Eycan be obtained by the difference discretization of equations (15) and (16) in the source region (Sun et al., 2013). We derive the difference equation of equation (17) in detail and, using the central difference, the derivative with respect to time is

By performing linear interpolation to Ez, we obtain

Substituting equations (18) and (19) into equation (17), we obtain the equation for Ezin the grid where the transmitting loop is located

The finite-difference spatial grids in this study are Cartesian grids, and the abovementioned transmitting current can be applied only on the gridlines. However, we need to rotate the minor transmitting loop in the roadways of coalmines. Therefore, it is necessary to study the loading of multidirectional minor transmitting loop into the Cartesian grids. The transmitted magnetic moment at any angle is M = IS, where I and S are the transmitting current and emission area, respectively. We assume that the components of M in the rectangular coordinate system are respectively Mx, My, and Mz, as shown in Figure 5, and by replacing the three components of M with the magnetic moments generated by the transmitting loop at three different directions, we obtain

where θ and φ are respectively the elevation and azimuth of the transmitting loop’s normal direction in the spherical coordinate system. Thus, the three currents in the transmitting loops are

The minor loop source for any transmitting direction can be calculated by using this method. Hence, the requirements for multidirectional detection in the roadways of coalmines by using minor loops are satisfied. In addition, we adopted the ramp step wave turn-off. The frequency spectrum of the ramp step wave is lower than that of the ideal step wave and thus it has larger wavelength. Therefore, it can better suppress the numerical dispersion.

Fig.5 Schematic of the multidirectional minor transmitting loop.

The following section is about the boundary conditions. The outer boundary of the full-space model consists of six faces. Unlike the half-space model, the transmitting source of the full-space model is located at the center of the model. If the grid range is sufficiently large, all the six boundary faces will be far from the transmitting source and the anomalies. In this case, the Dirichlet boundary conditions can be applied by setting the tangential electric field and normal magnetic field on all boundaries to zero. Because the air in the roadways is a nondissipative medium, the propagation in air of the electromagnetic waves satisfies the wave equation. If wedirectly use the wave equation to perform the iterations in air, we need to use a very small time step to satisfy the stability condition. The general equations for the iterations can be used as long as the resistivity of air is far higher than that of the other media in the model. Thus, not only can the general equations be used but also the computational precision can be guaranteed.

Verification of the calculation method

3D finite-difference modeling and analytical solutions in homogeneous full space

We used a multidirectional minor loop source based on the traditional finite-difference approach. This also improved the calculation of Bzin full space. In this section, we evaluate the applicability and precision of the method using a homogeneous fullspace model. Coincident loop devices were used. Both the transmitting and receiving loops were 2 m × 2 m square loops with normal directions parallel to the z-axis. The 1 A transmitting current was in the form of trapezoidal waves, with rising and falling edges of 1μs and pulse duration of 5 ms. A 3D nonuniform grid was used because nonuniform grids can reduce the computation time. Table 1 lists the details of the grids. The transmitting loop was located at the center of the model. In the rising and falling edges of the transmitting current, the time step was 1 ns and increased graduallywith time after turning off the current.

Fig.6 Comparison of the finite-difference solution and analytical solution for the homogeneous full-space model.

Table 1 Meshing parameters

Figure 6 shows the results for the finite-difference solution and the full-space analytical solution. The analytical solution is the wave response of the linear turn-off currentin homogeneous full space. To show the improvements in calculating Bz, Figure 6a shows the results obtained by the traditional method of calculating Bz. In Figure 6a, the curve of the electromotive force (EMF) in the finite-difference solution well coincides with that of the analytical solution, suggesting that the finite-difference modeling method used in this study is of high precision. Figure 6b shows the relative errors of the inductive potential produced by the proposed and the traditional method for calculating Bz, respectively. The maximum relative error of the traditional method is 25%, whereas that of the proposed method is below 5%.

Layered model

The layered model was used to verify the precision of the finite-difference method in this study. The model parameters are given in Table 2.

The transmitting loop was located in the middle of the third layer. Except for the geoelectrical model parameters, the rest were the same as in the homogeneous half-space model. Figure 7 shows the EMF decay for the homogeneous half-space model and the 1D digital filtering analytical solution. Initially, the decay of the finite-difference solution differs greatly from the analytical solution because the transmittingcurrent in the analytical solution is a negative step wave, whereas the current in the finite-difference solution is a ramp step wave. Over time, the decay curves of the finite-difference and the analytical solutions coincide.

Table 2 Model parameters

Fig.7 EMF comparison between the finite-difference solution and analytical solution for the layered model.

Full-space transient electromagnetic responses of different water-bearing structures

The precision of the computations was verified by using the full-space model in the previous section. The responses of water-bearing structures at variable distance and sizes were simulated and the resolution of the MTEM was studied.

Fig.8 Apparent resistivity of water-bearing structures at different distances vs. time.

Transient electromagnetic responses of waterbearing structure at different distances

A 3D low-resistivity body was set up ahead of the heading face to simulate the water-bearing structure. The earth model is shown in Figure 8a. The size of the transmitting loop, transmitting current, and transmitting waveform are the same as in the previous section. Figure 8a also shows the apparent resistivity of the water-bearing structure as a function of distance from the transmitting source. The apparent resistivity was calculated according to the method proposed by Yang et al. (2010). Owing to the effect of the turn-off time and the equation for the apparent resistivity, the initial apparent resistivity was high (Figure 8a). To better reflect the model conditions, the apparent resistivity data were corrected with the methods in the references (Raiche, 1984; Fitterman and Anderson, 1987). The corrected apparent resistivity is shown in Figure 8b. There are significant differences among the response curves of the water-bearing structures at variable distance from the transmitting source; moreover, the structure near the transmitting source shows a strong early response. When the waterbearing structure is located 30 m from the transmitting source, the anomaly in the curve occurs at 15 μs, whereas it occurs at 50 μs when the structure is 50 m from the source. In conclusion, as the distance of the water-bearing structure from the transmitting source increases, the weak low-resistivity anomaly appears later.

Transient electromagnetic responses of waterbearing structures of variable length

Figure 9 shows the apparent resistivity of waterbearing structures of variable length vs. the corrected turn-off time. From Figure 9, it can be seen that the large-volume water-bearing structures display more clearly the low-resistivity anomaly, which also lastslonger. The differences in the apparent resistivity of water-bearing structures of variable length suggest that the MTEM is a suitable detection method.

Fig.9 Apparent resistivity of water-bearing structures of variable length vs. time.

Layer thickness and the transient electromagnetic response

Figure 10 shows the transient electromagnetic response for layers of variable lengths as a function of the corrected turn-off time. It can be seen that the thickness of the layer where the transmitting source is located affects the apparent resistivity, especially in the beginning. The apparent resistivity increases with increasing layer thickness; however, the effect of the layer thickness decreases after 20 μs.

Fig.10 Apparent resistivity for variable length layers vs. time.

Full-space transient electromagnetic responses of coal-bearing strata

The transmitting and receiving devices in the MTEM are placed in the roadways of coalmines. Thus, the electromagnetic field distribution is affected by the electrical properties of the coal seam, and the received signals reflect the coal-bearing strata and the anomalies. Figure 11 shows the 3D model of the heading face in the mine. The roadway is located in the coal seam.

Fig.11 3D model of the coal seam and strata.

In the Huaibei mining area (Yu, 2007), the roof of the 6th main coal seam consists of sandy shale with resistivity set at 100 Ω·m. The coal seam has maximum resistivity of 400 Ω·m. The floor consists of siltstone with resistivity set at 150 Ω·m. The thickness of the coal seam was set at 10 m. The 4 m × 4 m roadway was located in the coal seam, and the resistivity of air in the roadway was set at 20000 Ω·m. Measurement points were placed at the heading face by using the sector detection method, as shown in Figure 1. Then, the 3D FDTD method was used to calculate the transient electromagnetic response at each measurement point. The size of the transmitting loop, transmitting current, and transmitting waveform were the same as in the previous section. The time–depth conversion method (Yu et al., 2008) was used to perform time–depth conversion of the data at each measurement point. Figure 12 shows the fan-shaped contour map of the calculated apparent resistivity.

In Figure 12, the horizontal axis represents the distance along the direction vertical to the axis of the roadway. The vertical axis represents the distance along the roadway direction and the origin of the coordinates corresponds to the position of the heading face in Figure 11. From Figure 12, it can be seen that the contour lines of the apparent resistivity are nearly circular with the heading face at the center and the apparent resistivity decreases with increasing radius. This can be explained as follows: because the roof and floor of coal seams have smaller resistivity than the coal seam, the apparent resistivity is mainly affected by the coal seam and thus the shallow parts have higher apparent resistivity. Because the diffusion range of the electromagnetic field into the surroundings increases with time, theeffect of the roof and floor on the received inductive signals increases and this causes the apparent resistivity to decrease with increasing detection range. From the forward numerical modeling results, one can draw the conclusion that the inherent differences in the electrical properties of the strata strongly affect the MTEM detection results. This helps to correctly distinguish the geological anomalies from the strata.

Full-space transient electromagnetic

Fig.12 Forward modeling of the coal-bearing strata.

responses of water-filled goaf

Transient electromagnetic responses ahead of the heading face

Water in goaves ponds and gathers in cavings and fractures formed by the deformation and fracturing of the overlying rock after the coal has been removed. The fractures in the overlying rock allow the water to move from the overlying aquifer to the cavings and fissures.

In coal mines, the water in goaves mainly exists in the cavings and its distribution and occurrence are difficult to control. There is no method yet to accurately calculate the height of the caving zones and water-flowing fractured zones. This is generally determined by using empirical equations. According to these equations, the height of the caving zones is generally 3–5 times greaterthan the mining height and the height of the fractured zones is generally 2–3 times larger than the height of the caving zones. Hence, the height of the caving zone formed by a 3 m mining height is about 12 m, and the height of the fractured zone is about 30 m. In this study, we assume that the water is mainly in the cavings and partly in the fractures; the lower parts of the fractured zone has better water-bearing capability. A 3D lowresistivity body was set up ahead of the heading face to simulate the water-filled goaf, as shown in Figure 11. Figure 13 shows the cross section on the y- and z-plane. Figure 14 shows the water distribution along the coal seam. The resistivity of the water-filled goaf was set at 0.5 Ω·m and the other parameters were the same as in the previous section. The detection was performed at the heading face toward the coal seam roof, the coal seam, and coal seam floor, as in Figure 15. In each direction, the measurement points were arranged using the sector detection method and the apparent resistivity sections inthe three directions were obtained by forward modeling using the 3D FDTD method, as shown in Figure 16.

Fig.15 Diagram of the detection directions.

Fig.13 Cross section on the yz-plane of the water-filled goaf.

Fig.14 Schematic of the water-filled goaf ahead of the heading face.

Figure 16 shows the apparent resistivity contour lines 30 m from the heading face. The contour lines are nearly circular with the heading face at the center, as in Figure 12. This is because the inductive signals received in the period corresponding to this depth range are not affected by the low-resistivity anomaly but by the strata. Forty meters from the heading face, the apparent resistivity is strongly affected by the low-resistivity anomaly. The low-resistivity anomaly is present in all three directions and is the lowest in the direction of the coalseam. From Figure 16b, it can be seen that the center of the lowresistivity anomaly coincides with that of the geological body. Because we use the 1D linear inversion method to calculate the apparent resistivity, there are differences

Fig.16 Apparent resistivity cross sections of the waterfilled goaf in front of the heading face: (a) toward the roof, (b) toward the coal seam, and (c) toward the floor.

Transient electromagnetic responses at the left side of the heading face

One major difference between the MTEM and the ground TEM is that the MTEM is affected by the full space. In other words, the MTEM detection results are affected by anomalies located along and opposite to the normal direction of the transmitting loop. Therefore, studying the full-space transient electromagnetic responses of anomalies in different directions is critical to interpreting the location of water-filled goaves. In between the shape of the anomaly and the geological body. Relative to the roof and floor, the area of the anomaly in the coalseam has the lowest resistivity and identical distribution to the model longitudinal distribution. The following conclusions can be drawn from Figure 16. The distribution of the apparent resistivity in all the three directions is jointly affected by the low-resistivity anomaly and the roof and floor of the coal seam. Near the heading face, the distribution of the apparent resistivity is mainly affected by the coal seamroof and floor, and the changes in apparent resistivity are reflected the resistivity difference between the two. Away from the heading face, the distribution of the apparent resistivity is mainly affected by the lowresistivity anomaly and not the geological strata.

Figure 17, the water in the goaf is moved to the left side of the roadway, whereas the other parameters are the same as in the previous section. Figure 18 shows the contours of the apparent resistivity in the direction of the coalseam.

Fig.17 Location of the water in the goaf water at the left side of the heading face.

As shown in Figure 18, there are two low-resistivity anomalies at the left and right side of the heading face. The anomaly area at the left side of the heading face is larger, has lower apparent resistivity, and is consistent with the position of the geological anomaly. However, the area of the anomaly at the right side is smaller, has higher apparent resistivity, and there is no geologicalanomaly in the area. The reason for the emergence of this false anomaly at the right side is the following. The MTEM is affected by the full space and the geological anomalies at both sides of the transmitting loop produce low-resistivity anomalies. When detection takes place at the right side of the roadway, the geological anomalies at the left side also produce anomalies. Thus, a false anomaly is detected at the right side. Because of the effect of the full-space to MTEM, the geological anomaly at the side of the roadway produces an anomaly at the other side. Clearly, the effect of the full space must be considered in the data interpretation to eliminate the interference of the false anomaly. This study lays the theoretical foundation for interpreting data from onesided detection of water in goaf.

Fig.18 Cross section of the apparent resistivity of the waterfilled goaf at the left side of the heading face.

Transient electromagnetic responses of waterfilled roadways ahead of the heading face

Many active coalmines are surrounded by abandoned mines within which are large numbers of abandoned roadways. The water gathered in these roadways is a potential safety hazard to the headings in coalmines. A water-filled roadway was placed in the coal seam of the model in Figure 11 and its distribution along the direction of the coal seam is shown in Figure 19. Its resistivity is set at 0.5 Ω·m and the other parameters were left unchanged.

The apparent resistivity of the water-filled roadway in Figure 19 is obtained in the direction of the coal bed and is shown in Figure 20. It can be seen from Figure 20 that there is a low-resistivity anomaly that corresponds to the location of the water-filled roadway, the anomaly curves as the water-filled roadway, and the location of the lowresistivity anomaly coincides with that of the water-filled roadway. In addition, there is a small low-resistivity anomaly at the left side. Based on the discussion in the previous section, this should be caused by the geological anomaly at the right side because of the full-space effect. Because 1D linear inversion method was used to calculate the apparent resistivity, there are differences between the shape of the anomaly and the geological body.

Fig.20 Cross section of the apparent resistivity of the waterfilled roadway ahead of the heading face.

Figures 16, 18, and 20 show the full-space transient electromagnetic apparent resistivity of the water in goaves at different locations and of different shapes. When using MTEM to detect the water in goaves, the location and range of the water ahead of the heading face can be determined to a first approximation by combining and comparing hydrogeological data, electromagnetic measurements, and forward modeling to help drilling.

Application

The results of the numerical modeling of the fullspace transient electromagnetic responses can be used in the MTEM detection of water in goaves. In this section, we use the method to detect water in goaf. To find out whether there are goaves and water-bearing structures ahead of the heading face, we used the MTEM to detect the water and ensure the safety of roadways. The coal seam belongs to the Lower Permian Shanxi Formation. The coal roof mainly consists of sandstone and mudstone. The floor of the coal seam mainly consists of sandstone. The coal seam is stable and dips between 4° and 6°. There are no faults in the workingface. The water mainly comes from fractures in the roof sandstone. Figure 21 shows the cross section of the apparent resistivity in the direction of the coal beds with the heading face interconnecting with the working face.

Near the heading face (Figure 21), the contour lines are affected by the roof and floor of the coal seam, and the contour lines of the apparent resistivity are nearly circular with the heading face at the center. The apparent resistivity contours decrease with increasing radius. In the range 50–90 m at the right front of the heading face, the apparent resistivity is less than 10 Ω·m and obviously lower than in the areas at the left side of the roadway and right ahead of the heading face, pointing to the obvious low-resistivity anomaly. The apparent resistivity in this low-resistivity anomaly is less than 10 Ω·m, whereas the apparent resistivity at the left side of the heading face and at the same distance from the heading face was more than 40 Ω·m. This suggests that the anomaly area at the right side shows higher reduction in resistivity and the rocksat the corresponding location have higher water-bearing capability. The geological data for the mine suggest that the coal seam is stable and there are no faults, collapse columns, and other geologic structures ahead of the heading face. Thus, it is speculated that the low-resistivity anomaly at the right side is due to the water in the goaf. By comparing the measurements with the modeling and detection results, it can be seen that the shapes of the contours are identical, which further confirms that the anomaly is due to the water in the goaf. The mine administration thus arranged for advanced drilling to test the low-resistivity area and detected water at 40 m from the heading face. The lowapparent resistivity anomaly detected by the MTEM was identical to water-filled goaf verified by drilling.

Fig.21 Cross section of the apparent resistivity obtained by using the MTEM.

The numerical modeling results were verified as well as the applicability of the MTEM indetecting waterfilled goaves. The experimental and numerical modeling results are identical and this further confirms that we can use a multidirectional minor loop source to perform advanced detection of water in goaves.

Conclusions

Based on the FDTD method, we considered the loading of minor transmitting loops in any direction into Cartesian grids by replacing the actual transmitting minor loop with multiple transmitting loops. We improved the method for calculating the z-component of the magnetic field by discretizing the electromagnetic field equation and its difference expression, and established the full-space 3D geoelectrical model. Then, using the sector detection method, we arranged the measurement points at the heading face, calculated the coal-bearing strata full-space transient electromagnetic responses, and performed forward modeling of waterfilled goaves of variable shapes at different locations.

The difference in the electrical properties between the coal seam and its roof and floor strongly affects the distribution of apparent resistivity, resulting in nearly circular distribution with the heading face at the center and the apparent resistivity contours decreasing with increasing radius. The low-resistivity spots in the apparent resistivity contours actually coincide with water-filled goaves.

When the goaf water was located in one side, there is a false low-resistivity small-scale anomaly on theother side owing to the full-space effect with increasing contours in this zone. In practice, it is necessary to combine the full-space effect and actual hydrogeological data to determine the location of water in goaves.

It was found that the experimental results were identical to the numerical modeling results. Therefore, the MTEM with multidirectional detection is an effective method.

Acknowledgments

We wish to thank Hu Wen-Bao, Hu Xiang-Yun, Cheng Jiu-Long, Shi Xian-Xin, and Cheng Jian-Yuan for their comments and suggestions.

Alford, R. M., Kelly, K. R., and Boore, D. M., 1974, Accuracy of finite-difference modeling of the acoustic wave equation: Geophysics, 39(6), 834-842.

Fitterman, D. V., and Anderson, W. L., 1987, Effect of transmitter turn-off time on transient soundings: Geoexploration, 24, 131–146.

Goldman, Y., Hubans, C., Nicoletis, S., et al., 1986, A finite-element solution for the transient electromagnetic response of an arbitrary two-dimensional resistivity distribution: Geophysics, 51(7), 1450-1461.

Jiang, Z. H., 2008, Study on the mechanism and technology of advanced detection with transient electromagnetic method for roadway drivage face: PhD Thesis, China University of Mining and Technology, Xuzhou.

Kaufman, A. A., and Keller, G. V., 1983, Frequency and transient soundings: Elsevier, New York.

Liu, Y., Wang, X. B., and Wang, B., 2013, Numerical modeling of the 2D time-domain transient electromagnetic secondary field of the line source of the current excitation: Applied Geophysics, 10(2), 134-144.

Oristaglio, M. L., and Hohmann, G. W., 1984, Diffusion of electromagnetic fields into a two-dimensional earth: A finite-difference approach: Geophysics, 49(7), 870-894.

Raiche, A. P., 1984, The effect of ramp function turn -off on the TEM response of layered earth: Exploration Geophysics, 15, 37- 41.

SanFilipo, W. A., and Hohmann, G. W., 1985, Integral equation solution for the transient electromagnetic response of a three-dimensional body in a conductive half-space: Geophysics, 50(5), 798-809.

Sun, H. F., Li, X., Li, S. C., et al., 2013, Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time: Chinese Journal of Geophysics (in Chinese), 56(3), 1049-1064.

Wang, T., and Hohmann, G. W., 1993, A finitedifference time-domain solution for three-dimensional electromagnetic modeling: Geophysics, 58(6), 797-809.

Yan, S., Chen, M. S., and Fu, J. M., 2002, Direct Timedomain numerical analysis of TEM fields: Chinese Journal of Geophysics (in Chinese), 45(2), 275-282.

Yang, H. Y., Deng, J. Z., Zhang, H., and Yue, J. H., 2010, Research on full-space apparent resistivity interpretation technique in mine transient electromagnetic method: Chinese Journal of Geophysics (in Chinese), 53(3), 651-656.

Yang, H. Y., and Yue, J. H., 2008, Response characteristics of the 3D Whole-Space TEM disturbed by roadway: Journal of Jilin University(Earth Science Edition)(in Chinese), 38(1), 129-134.

Yee, K. S., 1966, Numerical solution of initial boundary problems involving Maxwell’s equations in isotropic media: IEEE Transactions on Antennas and Propagation,, 14(3), 302-307.

Yu, J. C., 2007, Mine transient electromagnetic prospecting: Press of China University of Mining and Technology, Xuzhou.

Yu, J. C., Wang, Y. Z., Liu, J., and Zeng, X. B., 2008, Timedepth conversion of transient electromagnetic method used in coal mines: Journal of China University of Mining and Technology, 18(4), 546-550.

Zhdanov, M. S., Lee, S. K., and Yoshioka, K., 2006, Integral equation method for 3D modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity: Geophysics, 71(6), G333-G345.

Chang Jiang-Hao received his B.Sc. from China University of Geosciences (Wuhan) in 2012. Presently, he is a Ph.D. candidate in the China University of Mining and Technology, mainly working on the theory and applications of the mine transient electromagnetic method. E-mail: jhchang@126.com

Yu Jing-Cun is a professor at China University of Mining and Technology. His research interests are inmining geophysics and electromagnetic data processing.

Manuscript received by the Editor April 12, 2015; revised manuscript received September 4, 2016.

*This work was supported by the National Key Scientific Instrument and Equipment Development Project (No. 2011YQ03013307), the Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions, and Key Laboratory of Coal Resources Exploration and Comprehensive Utilization, Ministry of Land and Resources.

1. School of Resources and Geosciences, China University of Mining and Technology, Xuzhou 221116, China.

♦Corresponding author: Yu Jing-Cun (Email: yujcun@163.com)

© 2016 The Editorial Department of APPLIED GEOPHYSICS. All rights reserved.