Introduction

Front. Earth Sci.

Frontiers in Earth Science

Front. Earth Sci.

2296-6463

Frontiers Media S.A.

1778695

10.3389/feart.2026.1778695

Original Research

Mechanism of time-lapse electrical resistivity tomography (ERT) response to mining-induced fracture evolution in shallow coal seams: a coupled DEM-FEM study

Li et al.

10.3389/feart.2026.1778695

Yuteng

¹ ² * Writing - review and editing Methodology Writing - original draft Visualization Investigation Cheng

Jianyuan

² Validation Data curation Writing - original draft Wu

Zhengfei

² Validation Data curation Writing - original draft Zhao

Jiajia

² Data curation Validation Supervision Investigation Writing - review and editing

1 China Coal Research Institute, Beijing, China 2 CCTEG Xi’an Research Institute (Group) Co., Ltd., Xi’an, China

*Correspondence: Yuteng Li, chinasx_lyt@163.com

27 02 2026

2026

1778695

31 12 2025 21 01 2026 12 02 2026

2026

Li, Cheng, Wu and Zhao

https://creativecommons.org/licenses/by/4.0/

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.

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.

DEM-FEM coupling digital image processing (DIP) shallow coal seam time-lapse electrical resistivity tomography (ERT) water flowing fractured zone (WCFZ)

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.

section-at-acceptance

Geohazards and Georisks

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 ( k n ), shear stiffness ( k s ), and friction angle ( φ ). The constitutive relationship of the contacts follows the criteria below (Equation 1): Δ σ n = k n Δ u n Δ τ s = k s Δ u s I f τ s < τ max τ max = c + σ n ⁡ tan ⁡ φ I f t h e r e i s a s h e a r m o v e m e n t (1) where Δ σ n and Δ τ s are the normal and shear stress increments of the contact, respectively; Δ u n and Δ u s are the corresponding displacement increments; and c is the contact cohesion. When the normal stress of the contact exceeds the tensile strength ( σ t ), 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

Schematic diagram of UDEC meshing strategy. (a) Conventional Voronoi Tessellation (Uncontrollable Fracture Path); (b) Proposed “Brick” Meshing Strategy (Sedimentary Fabric and Anisotropy).

Diagram (a) shows a network of randomly sized and shaped polygons labeled "Random Polygon (Isotropic)," while diagram (b) depicts a staggered rectangular block structure with labeled vertical joints and horizontal bedding, highlighting tensile fracture and potential separation.

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

Flowchart of image processing and FEM mesh generation for mining-induced fracture networks. (a) Original image acquisition from the UDEC discrete fracture network; (b) Multi-phase segmentation and binarization based on Otsu’s thresholding; (c) Reconstruction of the heterogeneous FEM mesh.

Three-panel scientific graphic: (a) colored stratigraphic model with layered blocks showing geological deformation; (b) black-and-white binary image highlighting fractures, with inset histogram dividing matrix and fracture pixels by Otsu threshold; (c) finite element mesh with triangular grid overlaying blue fracture regions.

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 6000 × 1500 pixels for the 600 m × 150 m 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 I u , v , where u , v 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): G g r a y u , v = 0.299 R + 0.587 G + 0.114 B (2) where R , G , B 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 σ B 2 between the two media (fracture phase and matrix phase) (Equation 3): σ B 2 T = ω 0 T ω 1 T μ 0 T − μ 1 T 2 (3) where ω 0 and ω 1 are the probabilities of occurrence for fracture and matrix pixels, respectively, and μ 0 and μ 1 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 N row × N col and the physical resolution is dx × dy . 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): J = σ E (4) where σ is the conductivity of the medium (S/m), and σ = 1 / ρ is the resistivity ( Ω · m ).

Since the steady current field is an irrotational field, a scalar potential ϕ can be introduced such that E = − ∇ ϕ . Combined with the current continuity equation ∇ · J = − I δ r − r s , the Poisson’s equation for a point source field in a non-uniform medium can be derived (Equation 5): ∇ · σ ∇ ϕ = − I δ r − r s (5) where I is the current intensity (A), δ is the Dirac delta function, and r s 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): ∂ ϕ ∂ n = 0 (6) where n 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): ϕ r → ∞ = 0 (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

Schematic of the FEM computational domain, meshing, and boundary conditions.

Diagram showing a central study area with a fine mesh, surrounded by an infinite element domain labeled as IED. Surface boundary at the top is insulated, while side and bottom boundaries are marked as infinite elements with the infinity boundary at the bottom labeled V approaches zero.

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 ( ρ r / ρ w ≈ 100 ), 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

Geometric dimensions and boundary conditions of the numerical model.

Cross-sectional diagram showing layered geological strata: topsoil, overburden, main roof, immediate roof, coal seam, and floor. Dimensions labeled are length six hundred meters, height one hundred fifty meters, and coal seam depth one hundred twenty-three meters. Boundary conditions noted.

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 ( v y = 0 ), and the left and right boundaries are fixed in horizontal displacement ( v x = 0 ).

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/m³, and the densities of the roof and floor sandstones and siltstones ranged between 2400 and 2550 kg/m³.

TABLE 1

Physico-mechanical parameters of rock blocks.

Lithology	Density ρ (kg/m³)	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

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

Mechanical parameters of joints.

Stratum/Lithology	Normal stiffness K n (GPa/m)	Shear stiffness K s (GPa/m)	Friction angle ϕ j (°)	Cohesion C j (MPa)	Tensile strength σ t j (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

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 ρ w = 10 Ω · m , while the surrounding rock matrix and air were generalized as a high-resistivity background ( ρ r = 1000 Ω · m ). This setting ensures an electrical contrast of 1 0 2 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 1.0 × 1 0 − 5 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 T gray is set. The image matrix is traversed to classify pixels into two physical media, constructing an initial binary matrix B W 0 (Equation 8): B W 0 i , j = 0 , i f G r a y i , j < T g r a y 1 , i f G r a y i , j ≥ T g r a y (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

Comparison of the adaptive correction effect for the surface subsidence boundary. (a) Artifacts caused by global thresholding: the subsided air layer is misidentified as a low-resistivity fracture, leading to a non-physical “air short-circuit”; (b) Result of the surface contour recognition algorithm: the air-ground interface is correctly reconstructed, ensuring current flows physically through the subsurface strata.

Two-panel illustration comparing current flow modeling in subsurface environments. Panel (a) shows an incorrect approach with a red curved arrow labeled "Non-physical 'Air Short-circuit'" passing through a blue region above brown subsurface, indicating misidentified current flow in air. Panel (b) shows the correct approach with a green arrow labeled "Current Flow in Subsurface" passing through the brown layer beneath a dashed red line labeled "Insulating Air Layer," indicating current remains in the earth.

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., B W 0 i , j = 1 ) and define it as the true instantaneous surface elevation j surface at that position.

Air Layer Removal: Forcefully mark all pixel areas above the surface ( j < j 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 1

Limit 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 S w = 100 % ).

Assumption 2

Homogenized 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 1 0 2 ∼ 1 0 3 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: ρ w = 10 Ω · m ;

Rock Matrix/Air: ρ r = 1000 Ω · m .

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 A moves along the survey line, while the other current electrode B is located at infinity (or extremely far from the survey line). The potential electrodes M and N 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., T 0 , T 1 , . . . , T 20 ), 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 ( T 0 ), 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 ( T = 0 m) and three key nodes during the mining process ( T = 80 , 220 , 400 m) were selected. The evolution of overburden failure morphology is shown in Figure 6.

FIGURE 6

Evolution characteristics of overburden fracture networks at different mining stages. (a) First caving stage (40 m); (b) Separation development stage (120 m); (c) Fracture penetration to surface (220 m); (d) Goaf compaction and stability stage (400 m).

Four-panel illustration labeled a to d shows horizontal rock strata. Panel a depicts undisturbed layers. Panel b adds a fault and localized deformation. Panel c displays significant layer folding. Panel d exhibits extended folding and fault propagation.

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

True resistivity models of mining-induced rock mass reconstructed based on DIP. (a) Initial state (0 m); (b) First weighting stage (80 m); (c) Fracture penetration stage (220 m); (d) Stable mining stage (400 m).

Scientific figure with four panels labeled a, b, c, and d showing resistivity maps as color contour plots against height and distance. Panels reveal increasing complexity and development of blue contour features toward the center as the panels progress from a to d. Each plot uses a color bar ranging from blue, indicating low resistivity, to red, indicating high resistivity, with numerical values from zero to one thousand along the right side.

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 ( ρ r = 1000 Ω · m ). 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, ρ w = 10 Ω · m ). 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

Forward modeling results of apparent resistivity pseudosections during the mining process. (a) Initial background (0 m); (b) Weak response at first weighting (80 m); (c) Strong low-resistivity at penetration stage (220 m); (d) Strong on both sides and weak in the middle during the stable period (400 m).

Four sequential resistivity pseudo-section heatmaps labeled a through d show subsurface apparent resistivity variations, with depth on the y-axis and distance on the x-axis; a shows uniform background values, while b, c, and d depict increasing fracture-induced resistivity anomalies highlighted by linear yellow and green streaks, demonstrating the impact of fractures.

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 1000 Ω · m ), 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 x ≈ 50 m , depth z ≈ − 110 m ), a group of inclined low-resistivity bands extending toward the surface ( x ≈ 100 m ) appears. The decrease in apparent resistivity in this area is relatively small, with values dropping from the background of 1000 Ω · m to approximately 800 Ω · m . 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 ( x ≈ 50 ∼ 100 m ) deepens, and the apparent resistivity further decreases to approximately 600 Ω · m , indicating that the fractures there are fully penetrated and conductivity is enhanced. In the direction of working face advancement (bottom x ≈ 200 m extending toward surface x ≈ 300 m ), a new distinct inclined low-resistivity line appears, with a resistance value of approximately 700 Ω · m . 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 x ≈ 370 m to surface x ≈ 460 m area), a wide strong low-resistivity inclined anomaly band forms, with apparent resistivity significantly decreasing to 600 Ω · m . 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 900 Ω · m , 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 R is introduced. It is defined as follows: R x , z = ρ t x , z − ρ 0 x , z ρ 0 x , z × 100 % (9) where ρ t is the apparent resistivity at the mining time t , and ρ 0 is the background apparent resistivity before mining. The change rate imaging maps plotted based on this parameter are shown in Figure 9.

FIGURE 9

Imaging of apparent resistivity change rate during the mining process. (a) Local weak anomaly (80 m); (b) Expansion of strong damage (220 m); (c) “Strong-side, weak-center” differentiation and compaction recovery (400 m).

Three contour plots labeled a, b, and c display relative resistivity change rates versus distance and pseudo depth, each with a color scale from blue (negative) to red (positive); variation intensifies from plot a to c, with the strongest changes evident in plot c at greater distances and depths.

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 R ≈ − 20 % ∼ − 10 % . 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 R ≈ − 30 % ∼ − 40 % . 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 − 40 % (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 − 10 % or even approaches 0 % . 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 ( K n , K s ) primarily affects the magnitude of elastic deformation before failure. Sensitivity analysis indicated that varying stiffness within a reasonable range ( 1.0 ∼ 5.0 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” ( R ) as the imaging metric (Equation 9). Since R calculates the relative change before and after mining, the systematic errors caused by the absolute background resistivity ( ρ 0 ) 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 ( AM ) based on the target depth ( H ), typically following the rule of thumb AM ≈ H / 5 ∼ H / 3 . 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.

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.

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

Binley

Hubbard

S. S.

Huisman

J. A.

Revil

Robinson

D. A.

Singha

(2015). The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales. Water Resour. Res. 51, 3837–3866. 10.1002/2015WR017016

26900183

Cheng

J. L.

Peng

S. P.

(2014). Research progress and development trend of mining geophysical exploration technology. J. China Coal Soc. 39, 1742–1750.

Cheng

J. Y.

Zhu

M. B.

Wang

Y. 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.

Cheng

J. Y.

Wang

B. L.

Fan

(2022). Key technologies and typical application scenarios of coal mine geological transparency. Coal. Sci. Technol. 50, 1–12.

Cheng

J. Y.

Yan

Y. T.

(2025). Technical progress and direction of borehole geophysical exploration in coal mines. Coal Geol. Explor. 53, 1–21.

Comsol

A. B.

(2018). COMSOL multiphysics reference manual. Stockholm, Sweden: COMSOL AB.

Dahlin

Zhou

(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

Daily

Ramirez

Binley

LeBrecque

(2004). Electrical resistance tomography. Lead. Edge 23, 438–442. 10.1190/1.1729225

Dimech

Cheng

Chouteau

Chambers

Uhlemann

Wilkinson

(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

36285292

Fan

L. M.

(2017). Scientific connotation of water-preserved mining. J. China Coal Soc. 42, 27–35.

Huang

Q. X.

(2002). Ground pressure behavior and definition of shallow seams. Chin. J. Rock Mech. Eng. 21, 1174–1177.

Itasca Consulting Group (2014). UDEC (universal distinct element code) version 6.0 User’s manual. Minneapolis, MN, USA: Itasca Consulting Group.

Koestel

Kemna

Javaux

Binley

Vereecken

(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

Liu

Nie

Liu

Tian

Wang

(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

Wang

Ren

(2021). Mechanical properties and damage constitutive model of coal under the coupled hydro-mechanical effect. Rock Soil Mech. 42, 2.

Han

Zhao

(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

35568720

Y. T.

Cheng

J. Y.

J. J.

(2023). Advanced prediction method of mine DC resistivity method based on artificial neural network. Coal Geol. Explor. 51, 185–193.

Liu

Dai

Zhao

N. 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

J. J.

X. P.

Spitzer

(2010). Algebraic multigrid method for 3D DC resistivity modeling. Chin. J. Geophys. 53, 700–707.

J. J.

Wang

B. C.

Yan

(2019). Application progress of mine electrical method in coal seam mining failure and water hazard monitoring. Coal Sci. Technol. 47, 18–26.

J. J.

Wang

B. C.

D. S.

(2022). Application of mine resistivity monitoring system in water hazard prevention and control of coal mining face. Coal Geol. Explor. 50, 36–44.

Power

Gerhard

J. I.

Tsourlos

Soupios

Simyrdanis

Karaoulis

(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

Qian

M. G.

J. 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.

Qian

M. G.

Miao

X. X.

J. L.

(1996). Theoretical study of key stratum in ground control. J. China Coal Soc. 21, 225–230.

Tang

C. 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

Wang

Zhou

Zhang

Liang

(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

37024549

J. L.

Zhu

W. B.

Wang

X. 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.

Yue

Z. Q.

Chen

Tham

L. G.

(2003). Finite element modeling of geomaterials using digital image processing. Comput. Geotech. 30, 375–397. 10.1016/s0266-352x(03)00015-6

Zhang

Fan

Wang

(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

Zhao

Jin

Wang

(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.

Zheng

Wang

Lan

Ren

Tian

(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

39242957

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