Abstract
Introduction:
Dynamic monitoring of the Water Flowing Fractured Zone (WCFZ) is critical for preventing water hazards in shallow coal seams, yet mapping the complex spatiotemporal evolution of mining-induced fractures to electrical signals remains challenging.
Methods:
This study proposes a time-lapse electrical forward modeling strategy coupling Discrete Element Method (DEM) and Finite Element Method (FEM) via Digital Image Processing (DIP). A UDEC “Brick” meshing strategy was employed to simulate overburden mechanics, while a DIP-based pixel mapping technique reconstructed true resistivity models preserving geometric anisotropy for Pole-Dipole simulations in COMSOL.
Results:
Results reveal the “vertical initiation-penetration-compaction recovery” mechanism and its distinct electrical signatures. Specifically, the full penetration stage (220 m) forms short-circuit channels inducing strong low-resistivity anomalies. In the stable mining stage (400 m), the apparent resistivity section exhibits a typical “strong-side, weak-center” differentiation controlled by the “O-ring” theory. Quantitative imaging of time-lapse resistivity change rate confirms that boundary tensile zones maintain a high negative change rate of approximately -40%, while the central compacted zone recovers to about -10%.
Discussion:
This study validates the feasibility of using time-lapse electrical resistivity tomography (ERT) to quantitatively evaluate the “damage-recovery” state of the goaf, providing a theoretical basis for precise water hazard monitoring.
1 Introduction
Coal resources have long held a dominant position in China’s energy security strategy. With the strategic westward shift of China’s coal development, high-intensity mining of shallow coal seams in the western regions (e.g., Northern Shaanxi and Shendong areas) has become the norm. These mining areas are characterized by “shallow burial depth, thin bedrock, and thick loose layers” (Huang, 2002; Fan, 2017). Mining activities in such conditions easily disrupt the stress equilibrium of the overburden, causing the Water Flowing Fractured Zone (WCFZ) to develop rapidly upward and damaging groundwater resources. Unlike in deep mining, mining-induced fractures in shallow seams often penetrate directly through the overlying bedrock to surface loose aquifers, leading to severe roof water inrush and sand burst accidents (Xu et al., 2012). Therefore, elucidating the penetration evolution mechanism of mining-induced fractures in shallow seams and implementing real-time monitoring of associated water hazards are critical engineering requirements for ensuring safe production in western mining areas (Cheng et al., 2022).
Among various geophysical monitoring methods, the direct current (DC) resistivity method is widely used due to its high sensitivity to subsurface water bodies and ease of operation (Cheng et al., 2014; Li et al., 2015; Lu et al., 2019). Traditional static electrical methods are often constrained by geological background noise, making it difficult to identify weak mining-induced anomalies. In recent years, time-lapse electrical resistivity tomography (ERT) has demonstrated immense potential in fluid transport monitoring by capturing dynamic differences in the electrical properties of subsurface media at different times, effectively eliminating background formation effects (Daily et al., 2004; Binley et al., 2015; Dimech et al., 2022). For the specific geological conditions of western shallow coal seams, time-lapse ERT offers significant advantages in capturing the development channels of water-conducting fractures (Zhang et al., 2011; Zheng et al., 2024). However, the accurate interpretation of field monitoring data relies heavily on the precision of forward modeling. Clarifying the physical mapping mechanism between the spatiotemporal evolution laws of fracture networks and resistivity responses is the theoretical cornerstone for achieving high-precision inversion and early warning of water hazards.
Currently, numerical simulation studies on the resistivity response of mining-induced rock masses are mostly based on continuum mechanics frameworks (e.g., FEM or FDM). These methods typically employ equivalent medium theory, indirectly characterizing damage by weakening the physical property parameters of elements in the plastic zone. However, rock physics experiments indicate that the current conduction path in fractured media is mainly controlled by the connectivity, attitude, and geometric topology of micro-fractures, rather than a simple average of porosity. Traditional continuum models ignore the widespread discontinuity and anisotropy within the rock mass, failing to accurately depict the control of key water-conducting channels—such as the “Voussoir Beam” structure (Qian et al., 1996) and bed separation fractures—on the current field. This results in essential deviations between forward modeling results and real physical processes.
In contrast, the Discrete Element Method (DEM, e.g., UDEC) can explicitly simulate discontinuous mechanical behaviors such as bed separation, tensile fracture propagation, and block rotation. It has natural advantages in revealing the breakage of the “key stratum” in the overburden and the development mechanism of the “O-ring” (Itasca Consulting Group, 2014). Unfortunately, mainstream discrete element software lacks coupled electromagnetic field solvers, while mature electrical forward modeling software (e.g., COMSOL) struggles to handle the modeling of tens of thousands of dynamic random fractures. In recent years, with the development of transparent mine geology technology (Cheng et al., 2019), bridging the barrier between discrete mechanical models and continuous electric field models to achieve lossless transfer of fracture network geometric information and electrical response simulation has become an urgent interdisciplinary challenge to be solved.
Addressing these issues, this paper proposes a discrete-continuous coupling forward modeling strategy based on Digital Image Processing (DIP) (Yue et al., 2003). This method utilizes UDEC to refine the simulation of the dynamic failure process of the overburden in shallow coal seams, introduces image segmentation and pixel recognition technologies to extract high-precision geometric features of the fracture network, and maps them onto a finite element mesh to construct a non-uniform anisotropic electrical model. Subsequently, the Finite Element Method (FEM) is used to solve the full-process time-lapse resistivity response. This study aims to reveal the dynamic influence mechanism of “O-ring” evolution on apparent resistivity imaging from the perspective of “micro-fracture to macro-electric field” coupling, providing a solid physical basis for the refined monitoring of roof water hazards (Xu et al., 2012; Cheng et al., 2025).
2 Materials and methods
2.1 Discrete element mechanical modeling
The movement of overburden strata induced by coal mining is a complex, non-linear dynamic process involving stress redistribution, strata fracturing, rotation, and caving. Distinct from numerical methods based on continuum mechanics (such as FEM or FDM), sedimentary rocks in coal measures exhibit significant layered discontinuous characteristics. Their macroscopic mechanical behaviors and the formation of water-conducting channels are primarily controlled by discontinuities, such as bedding planes and joints (Tang, 1997). To explicitly simulate the dynamic evolution of mining-induced fractures, the Universal Distinct Element Code (UDEC) was selected as the mechanical simulation platform (Itasca Consulting Group, 2014).
2.1.1 Block-contact constitutive model
UDEC discretizes the rock mass into a collection of polygonal blocks cut by discontinuities. The interior of the blocks is divided into finite difference meshes to simulate the deformability of the medium, while normal and shear forces are transmitted between blocks via contacts. This method allows for large deformation, rotation, and even detachment of blocks, enabling the realistic reproduction of the dynamic evolution characteristics of the “three zones” in the roof under mining influence.
In the numerical calculations, rock blocks are treated as elastoplastic materials, using the Mohr-Coulomb yield criterion to describe their compressive-shear failure behavior (Liu et al., 2017; Li et al., 2021). The contacts between blocks employ the Coulomb-slip Joint Model, whose mechanical response is controlled by parameters such as normal stiffness (), shear stiffness (), and friction angle (). The constitutive relationship of the contacts follows the criteria below (Equation 1):where and are the normal and shear stress increments of the contact, respectively; and are the corresponding displacement increments; and is the contact cohesion. When the normal stress of the contact exceeds the tensile strength (), the contact opens, and both normal and shear stresses drop to zero, marking the formation of tensile fractures or bed separations.
2.1.2 “Brick” meshing strategy
Considering the significant layered sedimentary characteristics of coal measures, rock strata failure primarily manifests as interlayer bed separation and intra-layer fracturing. While the traditional Voronoi random tessellation method can simulate rock fragmentation, it fails to accurately reflect the regular fracture spacing and interlocking structures of sedimentary strata. Therefore, a deterministic “Brick” meshing strategy was adopted to construct the roof strata model.
This strategy configures horizontal bedding and vertical joints to intersect orthogonally, with the vertical joints arranged in a staggered pattern between layers. This configuration accurately captures the articulated structure of the “Voussoir Beam” and its instability and caving mechanism (Qian et al., 1996), As shown in Figure 1. Furthermore, the fracture network generated by this strategy inherently possesses geometric anisotropy, providing a high-precision geometric basis for the subsequent DIP-based anisotropic electrical modeling.
FIGURE 1
2.2 Mesh generation based on digital image processing
To achieve the lossless mapping of discrete element simulation results to continuous medium electrical models, this paper proposes a pixel-based finite element mesh reconstruction method (Comsol, 2018). This method can rapidly convert the geometric information of complex fracture networks obtained from UDEC simulations into heterogeneous electrical models required for finite element calculations. The specific flowchart is shown in Figure 2.
FIGURE 2
2.2.1 Image acquisition and preprocessing
First, the mining-induced fracture field calculated by UDEC is exported as a high-resolution RGB image. In this study, to ensure that the geometric details of the fracture network are preserved, the image resolution was set to pixels for the geological model. This corresponds to a physical spatial resolution of 0.1 m per pixel. This precision is sufficient to capture the macroscopic morphology of bed separations and vertical fractures while maintaining high computational efficiency for the subsequent Finite Element analysis. Let the original image be , where are the pixel coordinates. To eliminate color noise and reduce the computational dimension, the weighted average method is used to convert it into a grayscale image (Equation 2):where represent the brightness values of the red, green, and blue channels, respectively. The converted grayscale image preserves the geometric boundary information between the rock matrix and fractures, providing basic data for subsequent segmentation.
2.2.2 Multi-phase segmentation using Otsu’s method
Since the coal seam roof contains multiple media such as intact rock matrix, bed separation fractures, and water-conducting fractured zones, phase segmentation of the grayscale image is required. This paper employs the maximum between-class variance method (Otsu’s method) to automatically determine the segmentation threshold.
The algorithm traverses all possible gray levels to find the threshold that maximizes the between-class variance between the two media (fracture phase and matrix phase) (Equation 3):where and are the probabilities of occurrence for fracture and matrix pixels, respectively, and and are the corresponding mean gray values. Using this threshold, the image is discretized into a binary matrix, where 0 represents pores/fractures and 1 represents the rock matrix.
2.2.3 Pixel-to-element topological mapping
To construct the mesh model required for finite element calculations, this paper establishes a mapping relationship between the image pixel space and the finite element physical space. Each pixel in the image is treated as a four-node quadrilateral isoparametric element.
Node Coordinate Reconstruction: Assume the image resolution is and the physical resolution is . For any pixel, the coordinates of the four corner nodes in the physical space can be obtained through linear transformation.
Property Assignment and Model Export: Based on the phase labels in the binary matrix, the generated elements are grouped. The fracture element group is assigned low resistivity (e.g., resistivity of water), while the matrix element group is assigned high resistivity (e.g., resistivity of rock). Finally, a standard finite element mesh file containing Nodes, Elements, and physical properties is generated, which can be directly used for COMSOL electrical forward modeling (as shown in Figure 2).
2.3 Finite element resistivity forward modeling
2.3.1 Governing equations
In DC resistivity exploration, assuming the subsurface medium is isotropic and electrically stable, the electric field distribution follows the law of conservation of charge (Comsol, 2018). Under the quasi-static field approximation, the current density and electric field intensity satisfy Ohm’s law (Equation 4):where is the conductivity of the medium (S/m), and is the resistivity ().
Since the steady current field is an irrotational field, a scalar potential can be introduced such that . Combined with the current continuity equation , the Poisson’s equation for a point source field in a non-uniform medium can be derived (Equation 5):where is the current intensity (A), is the Dirac delta function, and is the position of the current source.
2.3.2 Boundary conditions
To obtain the definite solution of the above partial differential equation, reasonable boundary conditions must be established.
Surface Boundary: The ground surface is an insulating interface where no current flows out, satisfying the Neumann boundary condition (Equation 6):where is the normal vector of the boundary.
Infinity Boundary: At an infinite distance underground, the potential approaches zero, satisfying the Dirichlet boundary condition (Equation 7):
In finite element numerical simulation, to simulate the semi-infinite space and eliminate boundary effects caused by artificial truncation boundaries, an Infinite Element Domain was set up at the periphery of the model (Lu et al., 2010). This domain maps the finite geometric region to infinity through mapping functions, thereby precisely satisfying the zero-potential condition at infinity.
2.3.3 Implementation strategy
This study utilizes the AC/DC module of COMSOL Multiphysics software for solution (Comsol, 2018). The specific implementation process is as follows:
Mesh Import and Property Mapping: The generated finite element mesh file containing attribute information is directly imported into COMSOL. Through interpolation functions, the discrete binary resistivity attributes (fracture phase and matrix phase) are accurately mapped to each element in the computational domain, as shown in Figure 3.
FIGURE 3
Observation System Configuration: Regarding the observation system, the Pole-Dipole array was selected (Dahlin and Zhou, 2004). Compared to the Wenner array, this configuration offers greater depth of investigation and extremely high sensitivity to vertical plate-like low-resistivity bodies (such as water-conducting fractured zones), making it highly suitable for capturing the vertical development patterns of mining-induced fracture zones.
Physics Interface Settings: To simulate the actual Pole-Dipole observation setup, the current electrode is set as a point current source, and the potential at the infinite boundary is constrained to 0.
Solver Configuration: Considering the converted quadrilateral mesh and the high-contrast resistivity media (), a direct solver (such as PARDISO) is employed to solve the linear equation system, ensuring computational convergence and numerical stability.
3 Numerical simulation setup
3.1 Geological model
A two-dimensional discrete element numerical model was constructed based on the typical occurrence conditions of shallow coal seams in the Northern Shaanxi mining area. The geometric dimensions of the model are 600 m in length and 150 m in height, as shown in Figure 4. The simulated strata depth covers from the ground surface to the coal seam floor, with a coal burial depth of approximately 123 m, representing a typical shallow coal seam suitable for surface electrical detection.
FIGURE 4
The model is finely divided into 10 engineering geological layer groups from bottom to top. The coal seam is located at y = 20–27 m, with an average thickness of 7 m and a dip angle of 0°. The roof strata consist mainly of alternating layers of siltstone, fine sandstone, and medium-coarse sandstone, exhibiting typical sedimentary cycle characteristics. The topsoil layer extends directly to the top of the model, simulating the Quaternary loose layer.
To simulate the realistic semi-infinite space geological environment, the model boundary conditions are set as follows.
Displacement boundary. The bottom boundary of the model is fixed in vertical displacement (), and the left and right boundaries are fixed in horizontal displacement ().
Stress boundary. The top of the model is a free surface, with no additional equivalent load applied, directly simulating the surface topography.
Initial in situ stress. Considering the characteristics of shallow burial depth, the initial in situ stress field is dominated by gravity stress. The horizontal lateral pressure coefficient is set to 1.2, and the maximum horizontal principal stress is approximately 3–5 MPa.
3.2 Numerical implementation of “Voussoir Beam” structure and fracture evolution
The macroscopic failure mode of sedimentary rocks in coal measures is controlled by the spatial topological relationship between bedding and joints. After mining, the hard roof strata typically fracture to form neatly arranged rock blocks, which interlock above the goaf to form an articulated structure, namely, the classic “Voussoir Beam” structure. This mechanical process determines the connectivity and spatial distribution pattern of the water-conducting fracture network.
To realistically reproduce this key mechanical behavior in the numerical model and provide a medium model consistent with physical laws for subsequent electrical forward modeling, this paper adopts the aforementioned deterministic “Brick” meshing strategy.
3.2.1 Geometric topology reconstruction
In key strata such as the immediate roof and main roof, the traditional random Voronoi tessellation is abandoned in favor of a preset orthogonal cutting of horizontal bedding and vertical joints, with vertical joints distributed in a “staggered” pattern between layers.
Horizontal bedding: Simulates the bedding planes of sedimentary rocks, allowing for bed separation.
Vertical joints: Simulate the inherent weak planes of rocks, allowing for tensile fracturing.
This geometric configuration allows rock blocks to undergo free rotation and sliding after excavation unloading, thereby spontaneously forming an articulated arch structure with bearing capacity, effectively avoiding the issue of non-physical “loose body flow” often seen in conventional discrete element simulations.
3.2.2 Control of fracture anisotropy
The fracture network generated by this strategy possesses significant geometric anisotropy, which directly determines the response characteristics of the subsequent electrical model:
Horizontal direction: Prone to developing long-distance bed separation fractures, forming lateral conductive channels.
Vertical direction: Develops short-distance step-like broken fractures, forming vertical communication channels.
This specific geometric feature of fractures is the fundamental geological cause for the mining-induced rock mass exhibiting resistivity anisotropy response. Through the above strategy, the model not only ensures that the strata caving step in mechanical calculations matches field measurements but, more importantly, provides a high-precision geometric base map with clear physical significance for the DIP-based finite element electrical modeling.
3.3 Physico-mechanical parameters
The selection of parameters for the numerical model directly determines the reliability and physical authenticity of the simulation results. Therefore, the physico-mechanical parameters of the rock mass were determined with reference to relevant geological reports and laboratory test results from the typical shallow seams in the Northern Shaanxi mining area (Zhao et al., 2019).
However, rock masses in the field differ significantly from laboratory samples due to the presence of internal discontinuities. Considering the size effect of the rock mass, the intact rock strength parameters measured in the laboratory were reduced based on the Hoek-Brown criterion. Specifically, based on the Geological Strength Index (GSI) typical for the weathered bedrock and sedimentary strata in this region (estimated between 40 and 60), the macroscopic strength parameters were reduced to approximately 25%–30% of the laboratory values.
3.3.1 Block constitutive model and parameters
Rock blocks are treated as deformable bodies using the Mohr-Coulomb elastoplastic constitutive model. This model can well describe the yield and failure behavior of rocks under compressive-shear stress states.
To ensure scientific parameter management, rock strata were classified into four lithological categories based on similarity in lithological characteristics and mechanical properties (as shown in Table 1). The density of the coal seam was set to 1400 kg/m3, and the densities of the roof and floor sandstones and siltstones ranged between 2400 and 2550 kg/m3.
TABLE 1
| Lithology | Density ρ (kg/m3) | Bulk modulus K (GPa) | Shear modulus G (GPa) | Friction angle φ (°) | Cohesion C (MPa) | Tensile strength σt (MPa) |
|---|---|---|---|---|---|---|
| Siltstone/Mudstone | 2550 | 3.0 | 2.7 | 38 | 3.9 | 4.1 |
| Coarse sandstone | 2500 | 2.5 | 2.3 | 40 | 2.4 | 2.3 |
| Fine sandstone | 2400 | 2.4 | 2.1 | 36 | 2.0 | 1.7 |
| Coal seam | 1400 | 2.4 | 2.1 | 36 | 2.0 | 1.7 |
Physico-mechanical parameters of rock blocks.
3.3.2 Joint mechanical parameters
The macroscopic failure of the rock mass is mainly controlled by the mechanical properties of the joint surfaces. The contacts employ the Coulomb-slip Joint Model. To realistically reproduce the mining pressure characteristics, the joint parameters were calibrated against field observations. The assigned parameters (Table 2) ensure that the simulated periodic weighting step (approximately 15–20 m) and the surface subsidence factor are consistent with actual production data in the study area.
TABLE 2
| Stratum/Lithology | Normal stiffness (GPa/m) | Shear stiffness (GPa/m) | Friction angle (°) | Cohesion (MPa) | Tensile strength (MPa) |
|---|---|---|---|---|---|
| Main roof | 3.0 | 2.0 | 13 | 0.4 | 0.2 |
| Other sandstone/Siltstone | 2.6 | 1.4 | 15 | 0.1 | 0.05 |
| Coal seam | 2.0 | 1.0 | 18 | 0.3 | 0.2 |
| Immediate Roof1 | 1.4 | 0.7 | 18 | 0.3 | 0.2 |
| Floor | 2.6 | 1.4 | 15 | 0.1 | 0.05 |
Mechanical parameters of joints.
3.3.3 Electrical parameter selection
For the electrical model, the parameters were selected to represent the “Limit Water-Filling” scenario. Based on the hydrogeological conditions of the Northern Shaanxi coalfield, the resistivity of the freshwater-filled fractures was set to , while the surrounding rock matrix and air were generalized as a high-resistivity background (). This setting ensures an electrical contrast of orders of magnitude, which is typical for water inrush hazards in this region.
3.4 Mining procedure and electrical observation system
3.4.1 Step-by-step mining simulation strategy
To realistically reproduce the development process of the “two zones” in the shallow coal seam roof and the dynamic evolution of the overburden “Voussoir Beam” structure, the numerical simulation strictly follows the principle of “step-by-step excavation and progressive equilibrium”.
Mining Layout: The simulation adopts the retreating longwall caving mining method. The total length of the model along the strike is 600 m, of which 100 m–500 m is the mining face area. Protective coal pillars of 100 m are reserved on both the left and right sides of the model to eliminate the interference of boundary constraint effects on the mining influence zone.
Time Step Control: The mining step distance is set to 20 m, and the entire area is mined in 20 excavation steps. After each excavation command is executed, the program runs for 10,000 to 20,000 steps or until the unbalanced force ratio drops to to ensure that the caved rock mass reaches a quasi-static equilibrium state.
3.4.2 Electrical model reconstruction based on discrete fractures
To establish a quantitative mapping relationship between mining-induced rock mass damage and geoelectric field response, this paper applies the DIP mechanical-electrical coupling technology proposed in Section 2. The specific reconstruction process includes three key steps: binarization feature extraction of the fracture network, adaptive correction of the surface subsidence boundary, and physical mapping of electrical parameters.
3.4.2.1 Binarization feature extraction of fracture network
First, the displacement field and block distribution map calculated by UDEC are exported as high-resolution bitmaps. In the images, intact rock blocks appear as high gray values (or specific colors), while open bed separations and fractured cracks appear as low gray values (black or dark background).
Using image threshold segmentation technology, a gray threshold is set. The image matrix is traversed to classify pixels into two physical media, constructing an initial binary matrix (Equation 8):
At this point, the obtained matrix preliminarily characterizes the discontinuous structural features in the model but has not yet eliminated the air interference above the ground surface.
3.4.2.2 Adaptive correction of surface subsidence boundary
Since shallow coal seam mining causes significant non-uniform surface subsidence (such as bending subsidence basins and ground fissures), direct global threshold segmentation would misjudge the area originally belonging to “air” above the subsided surface as a fracture-like low-resistivity area (both are white in the binary map). In electrical forward modeling, if not distinguished, this would cause non-real short-circuit phenomena of current through the high-altitude air layer, severely masking the response of real underground fractures, as shown in Figure 5.
FIGURE 5
To ensure the physical authenticity of the model, this paper introduces a “surface contour recognition algorithm” to geometrically correct the binary model. The algorithm logic is as follows:
Column-wise Scanning: Scan each column of pixels in the matrix from top to bottom.
Boundary Locking: Identify the position of the first pixel in each column belonging to the “rock matrix” (i.e., ) and define it as the true instantaneous surface elevation at that position.
Air Layer Removal: Forcefully mark all pixel areas above the surface () as background media, retaining only the connected areas below the surface as effective fractures.
Through this step, geometric artifacts caused by large deformation are effectively eliminated, precisely defining the electrical computational domain of the “underground half-space”.
3.4.2.3 Electrical parameter mapping and model assumptions
Based on the corrected geometric model, definite resistivity physical properties are assigned to the media. To highlight the electrical contrast between mining-induced fracture water and surrounding rock media, and to facilitate subsequent change rate imaging analysis, the model is established based on the following two scientific assumptions:
Assumption 1Limit Water-Filling Assumption. Given the extremely shallow burial depth of the coal seam in the study area (123 m), it is assumed that the bed separation spaces and open fractures formed by mining are fully filled with water bodies (saturation ).
Assumption 2Homogenized Background Medium Assumption. Although there are certain differences in the resistivity of different lithologies (such as sandstone, mudstone, and coal seams) in actual strata, compared to the electrical contrast of orders of magnitude between water-filled fractures and surrounding rock, the electrical differences between strata are secondary factors. Moreover, the time-lapse resistivity change rate imaging method adopted in this paper can effectively eliminate the influence of static formation resistivity through background subtraction. Therefore, the model generalizes the undisturbed surrounding rock and surface air as a uniform high-resistivity background, setting only the water-filled fracture network as a low-resistivity anomaly body.Based on the above assumptions, the physical parameters are assigned as follows:Water-filled Fractures: ;Rock Matrix/Air: .A fixed-point monitoring mode is adopted. As the working face advances, full-section time-lapse resistivity observations are conducted on the surface (Power et al., 2015).
3.4.3 Surface time-lapse ERT system
To simulate the actual field detection process and obtain apparent resistivity responses in the numerical model, a high-density electrical observation system was established on the surface.
Survey Line Layout: The survey line is arranged at the top boundary of the model (surface), covering the entire strike length of the model. A total of 61 virtual electrodes are deployed, with an electrode spacing of 10 m. This layout meets the requirement for effective detection depth of the deep coal seam (burial depth 123 m) and ensures control accuracy for the lateral boundaries of the mining-induced fracture zone.
Observation Array: The Pole-Dipole Array is selected.
Array Principle: The current electrode moves along the survey line, while the other current electrode is located at infinity (or extremely far from the survey line). The potential electrodes and scan point by point along the survey line.
Selection Rationale: Compared to the Wenner array which is sensitive to horizontal layering, the Pole-Dipole array offers greater depth of investigation and extremely high sensitivity to vertical plate-like low-resistivity bodies (such as water-conducting fractured zones). This makes it capable of more clearly depicting the vertical development morphology of mining-induced fracture zones.
Time-Lapse Monitoring Strategy: A “fixed-point monitoring” mode is adopted. The surface electrode positions remain fixed. As the underground working face advances forward by 20 m (i.e., ), a full-section forward simulation data acquisition is performed on the surface. By calculating and comparing the apparent resistivity data at different times with the initial time (), the dynamic change characteristics of the underground electrical structure are captured, thereby achieving spatiotemporal imaging of mining-induced fracture evolution.
4 Results and analysis
4.1 Spatio-temporal evolution of roof fracture networks
To reveal the instability mechanism of the overburden “Voussoir Beam” structure and the development laws of the Water Flowing Fractured Zone (WCFZ) during shallow coal seam mining, a comparative analysis was conducted based on UDEC discrete element simulation. The initial state ( m) and three key nodes during the mining process ( m) were selected. The evolution of overburden failure morphology is shown in Figure 6.
FIGURE 6
Initial In-situ Stress Equilibrium State (Advancement 0 m). Before coal seam excavation, the model is in a state of in situ stress equilibrium. As shown in Figure 6a, the overburden bedding is clear, and interlayer contacts are tight, with no macroscopic tensile fractures or bed separation spaces developed. This state represents the original attributes of the geological body. Electrically, the dense rock skeleton appears as a continuous high-resistivity medium. This provides a clear “zero-damage” physical baseline for the subsequent identification of mining-induced fracture damage.
First Weighting and Formation of “Voussoir Beam” Structure (Advancement 80 m). When the working face advances to 80 m, the hanging area of the goaf reaches its limit, and the first weighting of the overburden occurs. The hard main roof fractures, forming a neatly arranged rock block structure. These fractured rock blocks interlock and articulate above the goaf, constituting a “Voussoir Beam” mechanical skeleton capable of bearing the load of the overlying strata. Accompanying the fracture of the main roof, obvious vertical broken fractures are generated within the rock strata. These fractures cut off the horizontal continuity of the strata, marking the transition of water-conducting channels from simple micro-bed separations to vertical penetration, providing an initial path for the subsequent short-circuiting of the electric field in the vertical direction (Qian et al., 1996).
Full Fracture Penetration and Formation of Surface Subsidence Basin (Advancement 220 m). As the working face advances to 220 m, overburden failure reaches its peak. As shown in Figure 6c, the WCFZ rapidly develops to its maximum height and triggers a significant surface response. The WCFZ is no longer confined to the interior of the rock strata but directly cuts through the top weathered rock and loose layers to reach the surface boundary. This marks the complete establishment of the “underground-surface” hydraulic connection, where surface water can easily rush into the mine along such penetrating fractures. The surface displays obvious step-like subsidence, and tensile cracks form at the edges of the goaf, destroying the integrity of the surface aquiclude (Wang et al., 2023; Li et al., 2022).
Goaf Compaction and Spatial Differentiation of Fracture Field (Advancement 400 m). When mining reaches the fully mined-out state (advancement 400 m), the overburden movement tends to stabilize, and the fracture field shows significant spatial differentiation, namely, the formation of the “O-ring” characteristic (Qian and Xu, 1998). In the central part of the goaf (within the range of x = 200–300 m in the figure), the caved and broken rock mass is gradually compacted under the immense self-weight of the overburden. Most of the originally open middle and low-level bed separation fractures close, and the contact tightness between rock blocks increases, leading to reduced water conductivity in this area. In contrast, on the open-off cut side (left) and the working face coal wall side (right), affected by the shear and cantilever tensile effects generated by the surrounding rock “abutment pressure”, the rock strata still maintain a high amount of dislocation, and fractures remain open. This mechanical state of “open fractures at both ends, closed fractures in the middle” forms a typical “strong-side, weak-center” structure.
4.2 DIP-based electrical model reconstruction
Using the aforementioned DIP technique, the discrete element calculation results of the four key mining stages in Figure 6 were mapped to continuous medium resistivity models, constructing a spatiotemporal evolution atlas of the true resistivity of the mining-induced rock mass, as shown in Figure 7. This figure intuitively reveals the complete evolution process of the underground electrical structure from “initial homogeneity” to “local development” and then to “full-line penetration and reconstruction” as the working face advances.
FIGURE 7
Initial High-Resistivity Background Field (Advancement 0 m). As shown in Figure 7a, before coal mining, the strata are in the in situ stress state, and no water-conducting fractures are developed. Based on the aforementioned “homogenized background medium” assumption, the model at this time appears as an overall uniform high-resistivity background (). This state constitutes the “zero-value” baseline for subsequent identification of mining-induced fractures and calculation of time-lapse change rates.
Initiation of Vertical Channels and Electrical Anisotropy (Advancement 80 m). As shown in Figure 7b, this stage corresponds to the first weighting of the main roof. The true resistivity model shows that with the formation of the “Voussoir Beam” structure, the originally continuous high-resistivity red background is cut by several low-resistivity bands (blue, ). The electrical structure at this time presents a significant step-like shape. Horizontally, interlayer bed separations provide lateral conductive paths; vertically, broken fractures communicate adjacent strata. This structure causes the macroscopic electrical properties of the rock mass to exhibit strong anisotropy, marking the transition of water-conducting channels from “interlayer incubation” to “vertical expansion.”
Characteristics of Fracture Penetration Stage (Advancement 220 m). As shown in Figure 7c (220 m), the electrical structure undergoes a qualitative change during this stage. The dark blue low-resistivity fracture zone is no longer confined to the interior of the rock strata but directly cuts through the top high-resistivity overburden to reach the surface, acting as an extremely efficient channel for current entry into the ground. This marks the formation of the optimal window period for electrical detection of water hazards.
Characteristics of Stable Mining Stage (Advancement 400 m). As shown in Figure 7d (400 m), in the fully mined-out stage, the electrical model presents typical “O-ring” closure characteristics. Both the open-off cut and working face sides maintain a dark blue low-resistivity state, corresponding to boundary tensile fractures. The color of the central area of the goaf changes from dark blue to red, indicating that with the closure of fractures, the resistivity of the rock mass rebounds significantly. This evolution process clearly validates the unique advantage of the DIP method in capturing the dynamic “damage-recovery” electrical response of the rock mass.
4.3 Apparent resistivity forward modeling response
The reconstructed non-uniform electrical model was imported into COMSOL Multiphysics software, and the Pole-Dipole observation system was adopted for forward calculation. By extracting potential data from surface nodes and applying geometric factor conversion, apparent resistivity pseudosections were plotted to analyze the apparent resistivity response laws of the mining-induced fracture zone, as shown in Figure 8.
FIGURE 8
Uniform Response Characteristics of Initial Background Field (Advancement 0 m). As shown in Figure 8a, before coal seam excavation, the strata are undisturbed by mining. Based on the homogenized background medium model established in this paper, the apparent resistivity section presents an overall uniform high-resistivity characteristic. The apparent resistivity values across the entire pseudosection are distributed evenly and remain stable (approximately ), appearing as a single red hue. There are no obvious vertical or horizontal resistivity gradient changes in the section, and the contour lines are sparse and gentle. This “pure” background field provides an ideal “zero-value” baseline for the subsequent precise identification of weak electrical anomalies caused by mining-induced fractures.
Weak Inclined Anomaly at First Weighting Stage (Advancement 80 m). As shown in Figure 8b, during the initial mining stage, no large-scale low-resistivity zone appears in the apparent resistivity section; only a weak electrical perturbation is observed on the open-off cut side. At the bottom left of the pseudosection (horizontal position , depth ), a group of inclined low-resistivity bands extending toward the surface () appears. The decrease in apparent resistivity in this area is relatively small, with values dropping from the background of to approximately . Analysis suggests that this corresponds to the initial breaking fractures at the open-off cut. Since the development scale of fractures is still small and the burial depth is relatively large (123 m), it only manifests as a weak anomaly in surface observations. The inclined morphology is a typical geometric distortion produced by the Pole-Dipole array when detecting vertical non-uniform bodies.
Bilateral Low-Resistivity Expansion at Penetration Stage (Advancement 220 m). As shown in Figure 8c, as the working face advances, the range of apparent resistivity anomalies expands significantly, exhibiting the characteristics of “new addition at the front, enhancement at the rear, and connection in the middle”. The color of the inclined anomaly band on the original open-off cut side () deepens, and the apparent resistivity further decreases to approximately , indicating that the fractures there are fully penetrated and conductivity is enhanced. In the direction of working face advancement (bottom extending toward surface ), a new distinct inclined low-resistivity line appears, with a resistance value of approximately . This corresponds to the currently moving “mining-induced fracture wavefront”. In the area between these two inclined low-resistivity lines, the overall apparent resistivity decreases, indicating that the overburden above the goaf is disturbed as a whole, and conductivity is generally enhanced.
“Two-Low, One-High” Differentiation at Stable Stage (Advancement 400 m). As shown in Figure 8d, in the fully mined-out stage, significant spatial differentiation occurs in the apparent resistivity section, and the anomaly morphology is finally finalized. The inclined low-resistivity band at the open-off cut on the left side remains. Simultaneously, at the new working face boundary on the right (bottom to surface area), a wide strong low-resistivity inclined anomaly band forms, with apparent resistivity significantly decreasing to . The most significant feature lies in the central area of the goaf. Unlike the overall decline in the 220 m stage, the apparent resistivity in this area shows a clear rebound, with values recovering to approximately , and some areas show little change. This distribution characteristic of “low resistivity at boundaries, recovery in the center” quantitatively reveals the compaction effect of the “O-ring” in the goaf. The inclined low-resistivity bands on both sides precisely lock the boundaries of the water-conducting fractures, while the resistivity rebound in the middle confirms the electrical recovery caused by fracture closure.
4.4 Imaging characteristics of time-lapse resistivity change rate
To eliminate the interference of formation background resistivity heterogeneity on mining anomaly identification and further highlight the spatiotemporal development characteristics of the water flowing fractured zone, the time-lapse resistivity change rate parameter is introduced. It is defined as follows:where is the apparent resistivity at the mining time , and is the background apparent resistivity before mining. The change rate imaging maps plotted based on this parameter are shown in Figure 9.
FIGURE 9
Local Weak Damage at First Weighting Stage (Advancement 80 m). As shown in Figure 9a, during the initial mining stage, the change rate imaging clearly captures the “germination period” damage near the open-off cut. A light blue inclined band appears in the lower left (open-off cut side), with a relatively small change rate amplitude of approximately . This quantitatively indicates that although the first weighting caused strata breakage, the fracture aperture and connectivity are not yet sufficient. The electrical structure of the rock mass is only slightly perturbed, and no large-scale conductive channel has formed yet.
Expansion of Strong Damage at Penetration Stage (Advancement 220 m). As shown in Figure 9b, as the water-conducting fractures penetrate the surface, the electrical damage of the rock mass reaches its peak. The blue anomaly area expands significantly, not only covering the open-off cut side but also extending forward with the advancement of the working face, forming a wide low-resistivity anomaly band. The color of the anomaly core turns deep blue, and the change rate reaches the minimum value of the entire process, approximately . This sharp drop in value intuitively reflects the full-line penetration of water-conducting fractures. At this time, the rock mass above the goaf is severely broken, and pore water filling leads to a drastic decrease in resistivity, indicating the moment of highest water hazard risk.
Damage Differentiation and Compaction Recovery at Stable Stage (Advancement 400 m). As shown in Figure 9c, in the fully mined-out stage, the change rate image exhibits a spatial differentiation characteristic of “strong-side, weak-center”, quantitatively revealing the self-healing effect of the rock mass. On both sides of the open-off cut (left) and the working face coal wall (right), the change rate remains at a high negative level of approximately (deep blue). This quantitatively proves that the tensile fractures under the action of the boundary cantilever beam remain permanently open, serving as strong water-conducting channels. In contrast to the 220 m stage, the color of the central area of the goaf fades significantly, and the absolute value of the change rate drops significantly to or even approaches . The data recovery from −40% to −10% strongly confirms the driving role of the rotation and compaction of the “Voussoir Beam” structure on fracture closure. This quantitative feature of “strong negative change at ends, weak negative change in the center” provides conclusive geophysical evidence for inverting the compaction state of the goaf using electrical methods (Li et al., 2023; Lu, et al., 2022; Koestel, et al., 2008).
4.5 Discussion on parameter sensitivity and robustness
Numerical simulations inevitably involve parameter uncertainties. To evaluate the reliability of the proposed approach, we analyzed the influence of key parameter variations on the electrical response.
4.5.1 Mechanical parameter uncertainty
The contact stiffness () primarily affects the magnitude of elastic deformation before failure. Sensitivity analysis indicated that varying stiffness within a reasonable range ( GPa/m) does not alter the macroscopic failure mode (i.e., the formation of the “Voussoir Beam” and “O-ring” structures), ensuring that the geometric basis for the electrical model remains valid.
4.5.2 Electrical parameter uncertainty
In actual engineering, the resistivity of mine water and rock strata may vary due to mineralization or lithology changes. However, this study utilizes the “Time-lapse Resistivity Change Rate” () as the imaging metric (Equation 9). Since calculates the relative change before and after mining, the systematic errors caused by the absolute background resistivity () are effectively cancelled out. Even if the absolute resistivity of the water body fluctuates, the characteristic pattern of “strong negative anomalies in water-conducting channels” remains distinct. Therefore, the “strong-side, weak-center” evolution law revealed in this study possesses high robustness and is not sensitive to minor perturbations in individual physical parameters.
4.6 Discussion on scale matching and field applicability
The successful application of the proposed DIP-ERT coupling strategy depends on the rational matching of three key scales: fracture scale, DIP pixel resolution, and electrode spacing.
4.6.1 Resolution matching mechanism
In this study, the DIP pixel size was set to 0.1 m, which serves as a “bridge” scale. It is larger than the microscopic aperture of a single crack (mm-scale) but significantly smaller than the spacing of the electrode array (10 m). This setting adopts an “Equivalent Porous Medium” approach at the pixel level—mapping dense micro-fracture zones into low-resistivity pixel clusters—thereby preserving the macroscopic anisotropy of the fracture network (e.g., the lateral continuity of bed separations) without requiring computationally prohibitive microscopic meshing.
4.6.2 Detection limits and scalability
The numerical results (Figure 8) indicate that the 10 m electrode spacing can effectively resolve the “Two-Zones” (caving zone and fractured zone) which typically span tens of meters vertically. However, for identifying thinner bed separations (e.g., <1 m), the current configuration may suffer from volume averaging effects.
4.6.3 Engineering implication
For field applications, the proposed workflow is scalable. Engineers should select the electrode spacing () based on the target depth (), typically following the rule of thumb . Meanwhile, the DIP forward modeling resolution should be kept at least 1/10th of the electrode spacing to accurately predict the theoretical electrical response and guide the inversion interpretation.
4.7 Validation and Engineering implications
To assess the credibility of the simulation results, we compared the identified electrical features with existing theories and field case studies.
4.7.1 Consistency with “O-ring” theory
The “strong-side, weak-center” apparent resistivity pattern identified in the stable mining stage (Figure 8d) is not an artifact but a geophysical manifestation of the mechanical “O-ring” theory (Qian and Xu, 1998). Classic rock mechanics studies indicate that mining-induced fractures form an “O-shape” circle, where boundary fractures remain open due to the cantilever effect, while the central goaf is re-compacted. Our simulated electrical response perfectly maps this mechanical process, validating that the DIP-based coupling strategy correctly captures the spatiotemporal evolution of the fracture field.
4.7.2 Rationality of the -40% threshold
The simulation suggests that a resistivity change rate of approximately −40% indicates a fully penetrated water-conducting channel. This value aligns with field observations in similar coalfields. For instance, time-lapse ERT monitoring in the Ordos Basin has reported that water inrush zones often exhibit resistivity decreases ranging from 30% to 60% relative to the background. Therefore, the −40% threshold derived in this study provides a reliable theoretical baseline for defining “High-Risk Zones” in early warning systems.
4.7.3 Engineering application strategy
In practical engineering, absolute resistivity values are often affected by geological noise. However, the relative change rate patterns revealed in this study are robust. We suggest that field engineers should focus on identifying the “U-shaped” or “O-shaped” high-magnitude anomaly bands at the working face boundaries, as these are the most probable pathways for roof water inrush.
5 Conclusion
This study proposed a time-lapse electrical forward modeling strategy coupling DEM and FEM to investigate the mechanical-electrical response mechanism of mining-induced rock masses. The main conclusions are as follows:
A pixel-based mesh reconstruction technique was developed to bridge the gap between discrete mechanics and continuous electrodynamics. By mapping the UDEC displacement field to a finite element mesh via digital image processing (DIP), this method effectively overcomes the limitations of isotropic equivalent continuum models, accurately capturing the geometric anisotropy and “finger-like” conductive channel characteristics of mining-induced fracture networks.
The simulation reproduced the “O-ring” evolution mechanism of roof fractures. The overburden failure process exhibits three distinct stages: vertical initiation during the first weighting (80 m), surface-underground short-circuiting during the fracture penetration period (220 m), and central compaction during the stable mining period (400 m). Mechanically, this forms a typical “tensile-boundary, compacted-center” structure.
The apparent resistivity response exhibits a characteristic “strong-side, weak-center” pattern. During the fracture penetration stage, a strong low-resistivity zone cuts through the entire overburden, indicating peak water inrush risk. In the stable stage, the anomaly differentiates: strong inclined low-resistivity bands persist at the boundaries (identifying the water-conducting fractured zone), while the resistivity in the central goaf recovers significantly.
Time-lapse resistivity change rate imaging enables quantitative evaluation of goaf compaction. The relative resistivity change rate recovers from approximately −40% (indicating strong damage and water-filled fractures at boundaries) to −10% (indicating compaction and closure in the center). This contrast provides a quantitative geophysical criterion for evaluating the “damage-recovery” state of the goaf.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
YL: Writing – review and editing, Methodology, Writing – original draft, Visualization, Investigation. JC: Validation, Data curation, Writing – original draft. ZW: Validation, Data curation, Writing – original draft. JZ: Data curation, Validation, Supervision, Investigation, Writing – review and editing.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This study was supported by the Natural Science Basic Research Program of Shaanxi Province (No. 2024JC-YBQN-0268; No. 2025JC-YBQN-410); the Science and Technology Innovation Fund of CCTEG Xi’an Research Institute (Group) Co., Ltd. (No. 2024XAYKF01; No. 2023XAYJS08). The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.
Conflict of interest
Authors YL, JC, ZW and JZ were employed by CCTEG Xi’an Research Institute (Group) Co., Ltd.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
BinleyA.HubbardS. S.HuismanJ. A.RevilA.RobinsonD. A.SinghaK.et al (2015). The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales. Water Resour. Res.51, 3837–3866. 10.1002/2015WR017016
2
ChengJ. L.LiF.PengS. P. (2014). Research progress and development trend of mining geophysical exploration technology. J. China Coal Soc.39, 1742–1750.
3
ChengJ. Y.ZhuM. B.WangY. H. (2019). Cascade construction of geological model of longwall panel for intelligent precision coal mining and its key technology. J. China Coal Soc.44, 2285–2295.
4
ChengJ. Y.WangB. L.FanT. (2022). Key technologies and typical application scenarios of coal mine geological transparency. Coal. Sci. Technol.50, 1–12.
5
ChengJ. Y.YanY.LiY. T. (2025). Technical progress and direction of borehole geophysical exploration in coal mines. Coal Geol. Explor. 53, 1–21.
6
ComsolA. B. (2018). COMSOL multiphysics reference manual. Stockholm, Sweden: COMSOL AB.
7
DahlinT.ZhouB. (2004). A numerical comparison of 2D resistivity imaging with 10 electrode arrays. Geophys. Prospect.52, 379–398. 10.1111/j.1365-2478.2004.00423.x
8
DailyW.RamirezA.BinleyA.LeBrecqueD. (2004). Electrical resistance tomography. Lead. Edge23, 438–442. 10.1190/1.1729225
9
DimechA.ChengL.ChouteauM.ChambersJ.UhlemannS.WilkinsonP.et al (2022). A review on applications of time-lapse electrical resistivity tomography over the last 30 years: perspectives for mining waste monitoring. Surv. Geophys.43, 1699–1759. 10.1007/s10712-022-09731-2
10
FanL. M. (2017). Scientific connotation of water-preserved mining. J. China Coal Soc.42, 27–35.
11
HuangQ. X. (2002). Ground pressure behavior and definition of shallow seams. Chin. J. Rock Mech. Eng.21, 1174–1177.
12
Itasca Consulting Group (2014). UDEC (universal distinct element code) version 6.0 User’s manual. Minneapolis, MN, USA: Itasca Consulting Group.
13
KoestelJ.KemnaA.JavauxM.BinleyA.VereeckenH. (2008). Quantitative imaging of solute transport in an unsaturated and undisturbed soil monolith with 3-D ERT and TDR. Water Resour. Res.44, W12411. 10.1029/2007wr006755
14
LiS.LiuB.NieL.LiuZ.TianM.WangS.et al (2015). Detecting and monitoring of water inrush in tunnels and coal mines using direct current resistivity method: a review. J. Rock Mech. Geotech. Eng.7, 469–478. 10.1016/j.jrmge.2015.06.004
15
LiB.WangZ.RenC. (2021). Mechanical properties and damage constitutive model of coal under the coupled hydro-mechanical effect. Rock Soil Mech.42, 2.
16
LiX.JiD.HanP.LiQ.ZhaoH.HeF. (2022). Study of water-conducting fractured zone development law and assessment method in longwall mining of shallow coal seam. Sci. Rep.12, 7994. 10.1038/s41598-022-12023-9
17
LiY. T.ChengJ. Y.LuJ. J. (2023). Advanced prediction method of mine DC resistivity method based on artificial neural network. Coal Geol. Explor. 51, 185–193.
18
LiuY.DaiF.ZhaoT.XuN. w. (2017). Numerical investigation of the dynamic properties of intermittent jointed rock models subjected to cyclic uniaxial compression. Rock Mech. Rock Eng.50, 89–112. 10.1007/s00603-016-1085-y
19
LuJ. J.WuX. P.SpitzerK. (2010). Algebraic multigrid method for 3D DC resistivity modeling. Chin. J. Geophys.53, 700–707.
20
LuJ. J.WangB. C.YanY. (2019). Application progress of mine electrical method in coal seam mining failure and water hazard monitoring. Coal Sci. Technol.47, 18–26.
21
LuJ. J.WangB. C.LiD. S. (2022). Application of mine resistivity monitoring system in water hazard prevention and control of coal mining face. Coal Geol. Explor. 50, 36–44.
22
PowerC.GerhardJ. I.TsourlosP.SoupiosP.SimyrdanisK.KaraoulisM. (2015). Improved time-lapse electrical resistivity tomography monitoring of dense non-aqueous phase liquids with surface-to-horizontal borehole arrays. J. Appl. Geophys.112, 1–13. 10.1016/j.jappgeo.2014.10.022
23
QianM. G.XuJ. L. (1998). Distribution characteristics of “O-shape” circle of mining-induced fractures and gas drainage technique in goaf. J. China Coal Soc.23, 466–469.
24
QianM. G.MiaoX. X.XuJ. L. (1996). Theoretical study of key stratum in ground control. J. China Coal Soc.21, 225–230.
25
TangC. A. (1997). Numerical simulation of progressive rock failure and associated seismicity. Int. J. Rock Mech. Min. Sci.34, 249–261. 10.1016/s0148-9062(96)00039-3
26
WangB.ZhouD.ZhangJ.LiangB. (2023). Research on the dynamic evolution law of fissures in shallow-buried and short-distance coal seam mining in lijiahao coal mine. Sci. Rep.13, 5625. 10.1038/s41598-023-32849-1
27
XuJ. L.ZhuW. B.WangX. Z. (2012). New method to predict the height of water flowing fractured zone based on position of key stratum. J. China Coal Soc.37, 762–769.
28
YueZ. Q.ChenS.ThamL. G. (2003). Finite element modeling of geomaterials using digital image processing. Comput. Geotech.30, 375–397. 10.1016/s0266-352x(03)00015-6
29
ZhangD.FanG.MaL.WangX. (2011). Aquifer protection during longwall mining of shallow coal seams: a case study in the shendong coalfield of China. Int. J. Coal Geol.86, 190–196. 10.1016/j.coal.2011.01.006
30
ZhaoC.JinD.WangH. (2019). Construction and application of overburden damage and aquifer water loss model in medium-deep buried coal seam mining in yushen mining area. J. China Coal Soc.07.
31
ZhengL.WangX.LanH.RenW.TianY.XuJ.et al (2024). Study of the development patterns of water-conducting fracture zones under karst aquifers and the mechanism of water inrush. Sci. Rep.14, 20790. 10.1038/s41598-024-71853-x
Summary
Keywords
DEM-FEM coupling, digital image processing (DIP), shallow coal seam, time-lapse electrical resistivity tomography (ERT), water flowing fractured zone (WCFZ)
Citation
Li Y, Cheng J, Wu Z and Zhao J (2026) Mechanism of time-lapse electrical resistivity tomography (ERT) response to mining-induced fracture evolution in shallow coal seams: a coupled DEM-FEM study. Front. Earth Sci. 14:1778695. doi: 10.3389/feart.2026.1778695
Received
31 December 2025
Revised
21 January 2026
Accepted
12 February 2026
Published
27 February 2026
Volume
14 - 2026
Edited by
Krzysztof Skrzypkowski, AGH University of Krakow, Poland
Reviewed by
Ling Zhu, Chengdu University of Technology, China
Chunwang Zhang, Taiyuan University of Technology, China
Updates
Copyright
© 2026 Li, Cheng, Wu and Zhao.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yuteng Li, chinasx_lyt@163.com
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.