Spatial and Temporal Distribution Pattern of Pre-Mining Grouting-Induced Microseismicity and Prediction of Water Inrush

1. Introduction

About 27% of China's proven coal resources are exposed to different degrees of water hazards, of which the floor water inrush hazard is particularly prominent in North China [1,2]. With the gradual increase in the depth of coal mining in China, the spatial relationship between the panel and the aquifer is more complicated, and the environment of geo-stress and permeability pressure is more variable, the severity and frequency of the occurrence of water inrush generally increase. Researchers have made many advances through statistical methods [3,4], mechanical analysis [5,6], numerical simulation [7,8] ,physical experiments [9,10], and other research methods, which have effectively revealed the mechanism of floor water inrush and proposed a variety of prevention and control methods. All of the above methods make too many assumptions, and the reality of the geological environment is very complex and variable, so there are often limitations for practical problems.

Microseismic monitoring is a non-destructive, real-time, dynamic, large-scale observation measurement technology [11,12], through monitoring the propagation of elastic waves in the stratum, it can reflect the performance of faults, joints and other discontinuous complex structures, water conduction channels and hydraulic source and other information [13,14]. With the rapid development of computer technology and three-dimensional visualization in recent years, it has been more widely used in the prevention of coal mine water inrush. Moment tensor inversion (MTI) theory was introduced in early seismic studies [15,16]. Microseismic information can be described mathematically using MTI theory, and so microseismicity due to rock fracture can be expressed as a linear combination of moments, which will be very helpful in understanding the type and mechanism of fracture [17].

In terms of application in the field, Jiang [18] obtained the measured parameters such as the activation pattern of the geological structure, the depth of the floor damage zone, the height of the fractured zone, and the size of the reasonable coal pillar through the three-dimensional display of the microseismic positioning results. Sun and Wang [19] established a prediction model of water-conducting channel based on microseismic monitoring signal mutation analysis by fitting the activation damage depth data from on-site microseismic monitoring. Jin [20] developed a method of identifying water-conducting channels based on the inversion of microseismic energy density and rupture connectivity of the rock layer, and verified it in the field by visual resistivity monitoring data. Zhou [21] utilized microseismic monitoring to classify the water-conducting channels into four types and find the weak areas of the grouting curtain, which well compensate for the weak links of water inrush. The above research results can be used to locate the position and types of water-conducting channels, classify the areas with high water-surge potential, and be applied in the prevention and control of water hazards on the floor and roof dynamic disasters in coal mines through microseismic monitoring. However, they all belong to the monitoring during the process of engineering disturbance and lack of timeliness, which means that violent water inrush accidents cannot be avoided if the water-conducting channel is connected to a water source with high water pressure and high water volume. There is a lack of research on the distribution pattern of microseismicity events induced by other man-made works (floor grouting, etc.) during the period when the working face has not yet been mined, and whether it can be used as an indicator for predicting water inrush. We refer to the microseismicity induced by grouting works prior to coal seam mining as “Pre-mining Grouting-induced Microseismicity, PMGIM”.

The purpose of this study is to propose a method that can predict the water inrush hazardous area in advance through the microseismicity-based information generated by the pre-mining grouting project. The mechanism of pre-mining microseismicity induced by grouting is explained by mechanical analysis and numerical simulation, and the temporal and spatial distribution characteristics and evolutional law of microseismicity induced by grouting area management in the working face during the pre-mining period are obtained by analyzing the measured data on site with the Wutongzhuang 182602 working face as the background. The prediction method of water inrush hazardous area in the working face is proposed to provide a theoretical basis for the detection of hydrological anomalies and potential water inrush hazardous area by using PMGIM.

2. PMGIM induction mechanism

2.1. Mechanical analysis

The primary objective of grouting treatment in water inrush-prone working faces is to utilize high-pressure slurry to displace water, thereby converting fractured aquifers into relatively impermeable barriers. The grouting treatment process can be categorized into three distinct stages: filling grouting, pressure-boost penetration grouting, and high-pressure expansion. During pressurized seepage grouting and high-pressure expansion grouting, the slurry initially permeates into fractures through seepage, gradually filling the fracture spaces [22]. As the amount of grouting increases and due to the incompressibility of the grout, internal pressure is exerted on the fracture walls, generating tangential tensile stress. When this tangential tensile stress surpasses the combined bonding force of secondary weak surfaces and the original tangential stress, the grout causes these weak surfaces to crack and further extends the fractures. The fractures continue to expand and extend until the grouting pressure can no longer induce further extension. Varying degrees of microseismic events are triggered by fracture initiation, propagation, and extension.

In this study, the rock mass is modeled as an elastic semi-infinite body. Assuming the rock mass to be a homogeneous and isotropic elastic medium, the stress state around the grouting hole can be simplified as a plane strain problem. Typically, grouting holes in these projects are circular; therefore, the established model of a flat circular hole is shown in Fig. 1. The rock mass surrounding the grouting hole experiences two principal stresses that are unequal and perpendicular to each other in the horizontal direction. Here, δ₁ represents the maximum principal stress, and δ₃ denotes the minimum principal stress. The rock mass around the circular hole is subjected to radial stress(δ_r), circumferential stress(δ_θ ), and shear stress (δ_rθ) at a distance (r) from the hole wall to the center of the grouting hole and at an angle (θ) with the X-axis.

Figure 1. Stress state of rock mass around grouting borehole.

Figure 1. Stress state of rock mass around grouting borehole.

According to elasticity theory, the stress state around an unlined hole in polar coordinates is given by：

(1)

Here, r is the radius of the drilling hole, δ is the radial distance, P₀ is the grout weight pressure, and P_d is the grouting pump pressure, which is zero during the filling and grouting stage.

When r =δ, the stress state around the hole is:

(2)

According to the maximum tensile stress failure criterion, the crack propagates along the path of least resistance. Specifically, when θ=0° or θ=180°, δ_θ = δ_θmin, and the crack extends in the direction perpendicular to the minimum principal stress. The mechanical criterion for the grouting pump pressure required for crack initiation is as follows:

(3)

Considering the increase in tangential stress due to radial seepage of slurry caused by the interaction between grouting pressure and pore water pressure P_i, the criterion for grouting pump pressure required for crack initiation is as follows:

(4)

(5)

Here, Δδ_θ represents the tangential stress increase, B represents the Biot constant of the rock mass, and μ represents the Poisson ratio of the rock mass.

Following crack initiation, as the grouting amount increases continuously and grouting pump pressure is continuously supplied, primary and secondary cracks in the rock mass expand further. The cracks induced by grouting are considered Griffith cracks, with an extension length of 2L [23].

When the fracture is horizontal, the stress intensity factor K_I-1 is

(6)

When the fracture is vertical, the stress intensity factor K_I-2 is：

(7)

When K_I exceeds the critical stress intensity factor K_IC, cracks continue to expand and propagate under the influence of pumping pressure, thereby triggering various degrees of microseismic events during this process.

(8)

2.2. Numerical simulation analysis of PMGIM

2.2.1. Numerical methods and modeling

Microseismic monitoring of the mine working face spans kilometers, covering a vast range. Attempting to simulate the entire floor numerically is impractical and could lead to significant errors or experimental failure due to overlooking small joints. Therefore, in this study, we selected a unit size of 160*160 cm after conducting several numerical simulation experiments, considering the actual grout diffusion radius and factors induced by micro-earthquakes during grouting. This size selection enhances the clarity of grouting engineering's fragmentation effects and rationalizes microseismic-intensive areas. The numerical simulation employs the parallel cementation model as its constitutive model, complemented by the smooth joint model for handling discontinuous structures like faults and folds.

Initially, a collection of rock particles is generated, and for scenarios involving discontinuous structures, smooth joints with a specified inclination azimuth are strategically inserted. Fluid exchange between domains is facilitated by pipe flow interacting with particles. Contact under the effective stress of particles results in compression or expansion, altering the permeability of the flow conduit. Fluid flow induces changes in pore pressure, thereby altering the effective stress on the particles and achieving fluid-structure coupling.

Initial confined water pressure (pore pressure) is incorporated into the model through boundary condition application. This is achieved using the servo control principle, where the sum of stiffness between the wall and the particle is computed, and then applied to the wall. Throughout the grouting process, servo control is consistently maintained, ensuring the stress boundary conditions remain unchanged. Initially, the water pressure transitions to free flow during the grouting process, followed by modification of grouting hole boundary conditions and application of a fixed grouting pressure of 24.24 MPa.

Model parameters are selected through a trial-and-error method [24], with specific details provided in Table 1 and Table 2. The grouting hole, located centrally, has a diameter of 152.4mm and contains a total of 17,405 particles, as depicted in Fig. 2

Figure 2. Model diagram.

Figure 2. Model diagram.

2.2.2. Influence analysis of primary fractures and structures on PMGIM

1.

Influence by primary fractures

Natural fractures are widely distributed in strata and are typically characterized using statistical methods. Fracture density serves as one of the key indicators used to describe these fractures. In two-dimensional models, fracture density (represented by ρ) denotes the number of intersecting joints per unit length. This section examines the impact of fracture development on grouting effectiveness. The physical and mechanical properties of the rock mass and ground stress levels remain consistent with those previously discussed. Crack size follows a power law distribution, and crack positions follow a uniform distribution. ρ is categorized into 3 levels: 6, 8, and 10, respectively.

Figure 3. Crack distribution characteristics under different crack density levels

Figure 3. Crack distribution characteristics under different crack density levels

Analysis of fluid injection simulation results reveals a strong correlation between grouting-induced crack distribution and primary crack density. Fig. 3 illustrates a significant increase in induced cracks with increasing joint density. At lower primary crack densities, induced cracks exhibit directional morphology, forming 3 main cracks surrounded by numerous smaller cracks, resulting in a 109.6% increase in total crack count. Increasing primary crack density alters crack morphology, increasing main cracks from 3 to 6 with significantly longer lengths, resulting in a 139.7% increase in total crack count. At a primary crack density of 10, crack numbers increase exponentially, with higher tensile stress concentrations at crack tips. Furthermore, many cracks merge with existing natural cracks, forming a complex crack network totaling 1878 cracks.

2.

Influence by geological structures

Grouting near geological structures results in distinct crack propagation patterns and varying microseismic events compared to more homogeneous areas. To address this issue, three experimental groups were devised to investigate how faults influence the grouting process at distances of 0.5, 0.6, and 0.7 m from the grouting holes, all under a consistent grouting pressure of 24.24 MPa. Fig. 4 presents the numerical simulation results.

Figure 4. Distribution characteristics of cracks at different distances from faults

Figure 4. Distribution characteristics of cracks at different distances from faults

Analysis of fluid injection simulation results reveals that the distribution of cracks induced by grouting correlates closely with the distance from the fault. As the distance from the fault decreases, the number of induced cracks significantly increases. At larger distances from the fault, induced cracks exhibit clear directional characteristics similar to the background scenario. Three main cracks are present, with their direction deflecting towards the fault. At a distance of 0.5m from the fault, the number of cracks increases dramatically to 1111, with concentration near the fault. The force chain diagram exhibits a pattern similar to the previous simulation model, with tensile stress concentration observed at the crack tips. However, the force chain on the fault appears nearly blank. The relatively weak surface of the fault hinders force transmission, causing grouting pressure to concentrate in the space between the fault and the grouting hole, resulting in thicker force chains and greater tension stress areas.

3. Spatio-temporal Distribution Characteristics of PMGIM

3.1. Overview of the Study Area

The 182602 working face is situated in the sixth mining area of Wutongzhuang Mine, with the coal seam buried at depths exceeding 950m. The strike length of the working face measures 950m, and it extends 256m along the dip. The working face predominantly extracts No. 2 coal, which features a stable and straightforward structure. The coal seam averages 3.2m in thickness and ranges in inclination from 8° to 23°. It is noteworthy that the Yeqing aquifer (YQ) beneath the coal seam exhibits reduced water abundance following grouting reinforcement. The Shanfuqing aquifer (SFQ) directly underlies the 2# coal seam, characterized by strong and heterogeneous water abundance. The Austrian ash aquifer (O₂) poses the primary challenge for water damage mitigation in the Wutongzhuang mining area due to its substantial water abundance, extensive distribution, considerable thickness, and robust connectivity.

According to drilling data from the grouting operations at the working face, grouting was applied above the shallow limestone aquifer and the Ordovician limestone aquifer, with a total grout volume of 7066.5 tons. The grouting holes were arranged in a "feather" or "belt" pattern, with adjacent branch holes spaced 50-60 meters apart. Details regarding the number of grouting holes and the specifics of the research area are illustrated in Fig. 5.

Figure 5. 182602 working face overview

Figure 5. 182602 working face overview

3.2. Analysis of microseismic correspondence in grouting disturbance area

Mining at the 182602 working face commenced on Nov. 7, 2018. The microseismic monitoring system was established as early as Nov. 2017 and has been continuously operational for 11 months, providing crucial data for pre-mining safety assessments. During this period, the regional grouting control project triggered varying intensities of microseismic events between the 2# coal seam floor and the Ordovician limestone karst aquifer, forming the focus of this study known as " PMGIM ".

Through detailed data analysis, we exclude microseismic interference induced by mining activities in the adjacent 182101 face and observe the distribution characteristics of microseismic events over different time periods:

In Fig. 6 (a), a total of 40 microseismic events were detected between January and February 2018, with 70% (28) concentrated within the grouting disturbance area.In subsequent periods (as depicted in Fig. 6(b) to (f)), spanning January through September 2018, the number of detected microseismic events rose steadily to 337. During this time, the proportion concentrated within the grouting disturbance area fluctuated slightly but generally stabilized between 62.0% and 67.9% as follows: 66.2% from January to February, 66.8% in April, 67.9% in May, 67.6% in July, and 62.0% in September.

Figure 6. Plane range of grouting disturbance zone

Figure 6. Plane range of grouting disturbance zone

Furthermore, microseismic events outside the grouting disturbance area exhibited scattered distribution without a clear clustering trend. The growth rate of these events decreased notably in the later stages of the grouting process (May 2018, September 2018), suggesting effective sealing of cracks as grouting advanced.

From the above analysis, it is evident that there is a strong correlation between primary microseismic events and the grouting disturbance area. Continuous seepage and fracturing caused by grouting lead to the destruction of water storage spaces and karst caves, disrupting the balance of the Austrian ash aquifer and releasing high strain energy. Simultaneously, this process promotes the expansion and propagation of primary and secondary cracks in the rock mass, triggering the clustering of microseismic events. Beyond the grouting disturbance area, as the influence on the Austrian limestone aquifer diminishes, the degree of fracture expansion decreases, resulting in a discontinuous microseismic distribution.

3.3. PMGIM distribution law of 182602 working face

The occurrence of microseismic events indicates that the stress of the source rock mass decreases and the elastic deformation energy is released within a certain time and space. In order to facilitate the analysis of the time distribution, spatial distribution and energy characteristics of PMGIM on the working face of 182602, combined with the stop-mining time of the adjacent working face 182101 (February 22, 2018), The 182602 face PMGIM can be divided into two periods.

Period 1: involving January and February 2018, affected by the mining of the adjacent 182101 working face and the regional grouting treatment project.

Period 2: involving March, June and September 2018, affected by regional grouting management projects.

3.3.1. PMGIM Plane Distribution

The kernel density estimation method centers on the location of specific element points and reflects the clustering state of elements in the current space by calculating the density value of each grid under varying search radii. This method is suitable for analyzing the planar distribution patterns of PMGIM, and the calculation formula is as follows [25]:

(9)

Here, P_n(x) represents the kernel density value of the x-th microseismic event in the study area; R is the search radius; n is the number of all microseismic events within the search radius; K_i is the weighting function; (x-x_i)²+(y-y_i)² is the square of the distance between two microseismic events at (x,y) and (x_i,y_i).

Figure 7. Nuclear density distribution of PMGIM events at 182602 working face

Figure 7. Nuclear density distribution of PMGIM events at 182602 working face

Fig. 7 illustrates the kernel density cloud map of microseismic events at the 182602 working face. During period 1, the core density of PMGIM events was centered on the Zhonglicun syncline axis, forming a "V"-shaped, two-center continuous pattern. The high-frequency core density zcenters "A" and "B" were positioned near the main grouting hole and the central region of the 182101 working face, respectively. As a result of continuous grouting, the range of core density center "A" in February 2018 expanded significantly. In period 2, the planar distribution of core densities for microseismic events exhibited phased changes. In March 2018, the nuclear density values displayed a two-center discontinuous pattern extending from the northwest to the southeast. By June 2018, the nuclear density values revealed a central discontinuous pattern oriented from the northwest to the southeast. Compared to March 2018, the density value of core "A" rose by 19.2%, while the density value of core "B" declined by 41.3%. By September 2018, the nuclear density values exhibited a discontinuous layout from northwest to southeast, with the two centers disappearing. Compared to June 2018, the core density value of center "A" decreased by 45.1%, and that of center "B" decreased by 41%. During period 2, there was virtually no microseismic event response in the 182602 working face, indicating effective sealing of the fractures.

3.3.2. Frequency - energy characteristics of PMGIM

A total of 900 PMGIM events were recorded on the floor of the coal seam between June 10 and August 20, prior to mining at the 182602 wroking face. Fig. 8 illustrates how these events are categorized into four distinct stages based on their frequency and energy characteristics, as follows:

Phase I: characterized by low frequency and low energy, with a daily frequency <15, total energy <10⁵J, and single event energy <10³J.

Stage II: characterized by high frequency and low energy, with a frequency ≥15, total energy ≤2×10⁵J, and 95.8% of events with energy <10³J.

Phase III: characterized by high frequency and high energy, with a frequency ≥20, total energy >5×10⁵J, and a short duration.

Stage IV: transition stage with large frequency fluctuations, small total energy, maximum <2×10⁵J, and 87.4% of events with energy <10³J.

Since July 10, the microseismic events have repeated the aforementioned phase cycle.

Figure 8. Characteristics of daily frequency-energy variation of the floor microseismic events

Figure 8. Characteristics of daily frequency-energy variation of the floor microseismic events

3.3.3 Characteristics of agglomeration law

Prior to mining in the working face, PMGIM dynamics adhere to three major patterns: sporadic, clustered, and staged. Sporadic microseismic activity suggests non-uniform distribution closely related to regional stress variations. For instance, the initiation and cessation of grouting activities significantly impact microseismic frequency, facilitating real-time monitoring of grouting effects. Clustered microseismic events are concentrated in geologically fragile areas, such as fracture zones, indicating tectonic stress anomalies and potential water conductivity risks. The staged pattern of microseismic energy release progresses through three phases: "low frequency low energy - high frequency low energy - high frequency high energy," consistent with the grouting process. This pattern aids in predicting dangerous periods in the working face, particularly during the "high frequency high energy" stage, when vigilance for water gushing risks is imperative.

4. Prediction of water inrush risk based on PMGIM

4.1. Agglomeration evolution of PMGIM

The more developed the natural fractures in the rock mass (comprising various defect structures), the more easily microseismic events are triggered by grouting, leading to the formation of abnormal PMGIM accumulation areas. The source distribution and concentration degree of PMGIM events serve as crucial discrimination indices. Therefore, the Moran index is utilized to measure the spatial clustering and concentration degree of microseismic event attributes within the working face, thereby identifying potential water inrush areas [26,27].

4.1.1. PMGIM agglomeration mechanism

The variation trend of PMGIM is primarily influenced by structural polarization effects and the diffusion effect of grouting. The structural polarization effect, such as faults at the 182602 working face and the Zhonglicun synclinal axis, indicates that the injected slurry increases pore pressure in the formation. Consequently, the corresponding Mohr circle of the rock mass shifts leftward, moving the Mohr-Coulomb circle below the envelope and bringing the Mohr circle closer to the strength line. When the shear stress applied to section MN exceeds the shear strength of the rock mass (τ_s), shear failure occurs, leading to fracture expansion and triggering microseismic events (see Fig. 9).

The more developed the natural structures are, the more easily microseismic events are stimulated during the grouting process, forming the HH region characterized by high microseismic activity and concentration. The diffusion effect of grout implies that as the grouting process continues, cracks in numerous areas are effectively sealed and rock layers are strengthened. This makes it increasingly difficult to excite microseismic events, resulting in a significant reduction in their occurrence.

Figure 9. Diagram of pore pressure and effective stress of rock mass

Figure 9. Diagram of pore pressure and effective stress of rock mass

4.1.2. Moreland index analysis of PMGIM in 182602 working face

The analysis route for identifying the dangerous water inrush areas on the working face using the Moran index in this study is illustrated in Fig. 10.

Figure 10. Analysis of microseismic events based on Moreland index.

Figure 10. Analysis of microseismic events based on Moreland index.

The Moreland index value of PMGIM event frequency at the 182602 working face is 0.59, indicating a positive spatial correlation and a clustering pattern with 99% confidence. The P value of PMGIM events is 0.006, and the Z value is 11.64, indicating the presence of extreme values (statistically significant hot spots) within the working face. This suggests that PMGIM events are highly locally clustered. Since PMGIM events are concentrated in areas where natural fractures of the rock mass are more developed, these areas with higher concentrations are potential water inrush danger zones.

4.2. Hazard prediction for the 182602 working face

Firstly, the 182602 working face is meshed and coded, and a total of 90 units are set. The local Moreland index was calculated to determine the specific location of the anomaly area of PMGIM events, and the frequency value of PMGIM events was assigned to each unit, as shown in Fig. 11.

In the initial phase of the working face grouting project (January 2018), a higher concentration of HH units was observed, clustered in the 400-800m section and adjacent to the return air roadway densely. PMGIM events were predominantly distributed across the working face (with minimal LL units), indicating a poor grouting effect and a significant risk of water inrush.

By the intermediate phase of the working face grouting project (April 2018), HH units remained densely distributed in the 400-800m section near the return air roadway, but their numbers had significantly decreased. HL and LH elements appeared randomly, while the discrete elements of PMGIM events gradually decreased. This indicates that the grouting project successfully sealed some of the cracks in the working face, improving the overall grouting effect.

In the final phase of the working face area treatment project (October 2018), HH units significantly decreased, although their aggregation location remained unchanged. Basic PMGIM event distribution (LL units) covered most areas of both the inner and outer sections of the working face. At this stage, LL units without microseismic events largely covered the working face, though some HH units persisted, indicating that while grouting reinforcement was largely complete, there remained a risk of water inrush in the HH regions.

According to the PMGIM distribution results using the local Moran index, regions at high risk of water inrush were identified. One area includes units 81, 82, 83, 64, and 65HH, located 400-550m along the return air lane (near the F602-4 fault area). The other area includes units 85, 86, 87, 68, 69, 70, and 51, located near the Zhong Lishun syncline axis, as illustrated in Figure 11.

4.3. Electrical survey

The network parallel electrical method is founded on analyzing and interpreting apparent resistivity data of rock mass, where areas of low resistivity typically denote water-rich anomalies or water-bearing structures. Fig. 12 depicts the results of parallel electrical geophysical exploration on the 182602 working face. Along the strike, four anomalous areas (AA1, AA2, AA3and AA4) were identified in the outer section of the working face, confirmed later by drilling to have significant water presence, closely matching the high-risk zones for severe water inrush identified by Moran index analysis. Abnormal microseismic events persisted in area AA4 during mining operations, escalating the risk of water inrush. Subsequently, a drainage borehole was constructed, discharging water at a maximum rate of 130m³/h to prevent water inrush. This underscores the primary PMGIM event as pivotal in delineating areas prone to water inrush.

5. Conclusion

In this paper, the concept of Pre-mining Grouting-induced Microseismicity was innovatively proposed. Mechanism of PMGIM is revealed through mechanical analysis; subsequently, the influence of multiple factors on PMGIM is deeply explored by numerical simulation. Taking the 182602 working face in Wutongzhuang coal mine as an example, the temporal and spatial distribution characteristics of PMGIM in this area were summarized in detail. On this basis, we introduced the Moran index to predict the potential location of water-rushes danger in the working face. The main conclusions are as follows:

The interaction between the high-pressure slurry and the surrounding rock body causes tangential tensile stresses on the inner wall of the cracks, which lead to crack initiation and expansion of the cracks. Numerical simulation results show that discontinuous geological structures, such as joints and faults, degrade the mechanical properties of the formation and change the crack initiation and expansion pattern caused by the grouting itself, which is the main reason for the surge of microseismic events.

There is a good correspondence between the PMGIM event and the grouting disturbed area, and the distribution of planar kernel density has an obvious stage, and the stage pattern of energy release of the PMGIM event is consistent with the grouting process, which has gone through the three phases of “low-frequency low-energy, high-frequency low-energy, and high-frequency high-energy”. This spatial and temporal distribution characteristic is helpful to predict the danger period of water breakout in the working face, especially in the stage of “high-frequency and high-energy”.

The location of natural fracture development is more likely to be affected by grouting action, thus stimulating the occurrence of PMGIM events and the formation of PMGIM clustering areas. The Moran index analysis was introduced to analyze the evolution law of four types (HH-type, HL-type, LH-type, LL-type) of PMGIM event aggregation zones, and it was determined that the area with higher degree of PMGIM aggregation is the potential water-surge area.