#### Abstract

In order to extract gas accurately in Yue Nan coal mine and prevent gas over limit and gas accidents, a combination of theoretical analysis and numerical simulation was used to investigate the height of the overburden caving zone and fracture zone in the working face, using 215102 working face as the engineering background. The results show three key strata in the 215102 working face, namely siltstone, sandy mudstone, and sandy mudstone. The empirical formula’s calculation results are in good accord with the findings of the theoretical analysis, which indicates that the height of the water-conducting fracture zone created by mining in 215102’s working face is 87.35 m. The plastic zone, stress distribution and displacement variation of the model overburden of the working face were analyzed separately in the numerical simulation. The results of the plastic zone simulation show that the caving zone’s maximum height is about 15.96 m, and the fracture zone’s maximum height is approximately 82.39 m. The stress distribution shows that the caving zone’s greatest height is around 13.65 m, while the fracture zone’s maximum height is roughly 77.24 m. The amount of overburden subsidence proves that the caving zone’s maximum height of about 10.08 m, and the fracture zone’s maximum height of approximately 82.69 m. The height of the overburden caving zone and fracture zone of 215102 working face are ultimately found to be 14.02 m and 76.42 m, respectively, based on theoretical analysis, empirical formula calculation, and numerical simulation findings.

#### 1. Introduction

When the coal seam is not mined, the working face and its surrounding rock seam are in stress equilibrium, but as coal seam mining continued, the overburden’s stress balance was disturbed, causing the overburden to break and collapse and the destructive stress to rebalance. The continuous mining activities in the working face lead to an abundance of fractures in the upper rock stratum. The “three overburden zones”—the caving zone, the fracture zone, and the continuous deformation zone—are created during coal mining when the overburden strata collapse and fracture. The water-conducting fracture zone (WCFZ) also refers to the caving zone and the fracture zone [1–3]. For high gas mines, the WCFZ in the overburden has become the main area of gas migration and desorption [1]. Gas drainage in this region is a crucial tool for achieving okay coal mining and reducing financial losses. Therefore, studying and determining the height of the overburden WCFZ can effectively prevent accidents from occurring in the mine and guide safe mine production. Different nations compute the height of WCFZ using various techniques. For example, the British Coal Authority issued the regulations on underwater coal mining in 1968, which stipulated the composition of overburden, thickness, coal seam mining thickness, and mining methods [4]. The former Soviet Union published the method guide to determine the height of WCFZ in 1973, and according to the depth of the soil layer, how frequently mining happens, and other considerations, the relevant coal mining underwater body requirements released in 1981 explain the safe and appropriate mining depth [5]. In China, scholars usually specify the height of WCFZ in two ways, a determination technique due to the placement of important strata in the overburden and an empirical calculation based on years of statistical data [6, 7].

Numerous academics have studied the height of development and the law of the WCFZ for a very long period, and their work has produced numerous findings that have received widespread acclaim in the field. The hypothesis of “key strata” was proposed by academics to Qian et al. [8] for the hard rock strata that have a significant controlling function in rock strata activities. The position of the key strata may be used to predict the height of the WCFZ, and Xu et al. [7, 9] established a practical method to identify the key strata by this theory. Based on the theoretical calculation of the height of the WCFZ in the key strata, Wang et al., Lu et al., and Cao et al. [10–12] studied whether there are fractures, fracture height, and formation settlement deformation with different fracture penetration rates in the continuous deformation zone of the longwall working face. In the overburden of the Donetsk coalfield in Ukraine, Palchik [13–15] used the drilling peep method to examine the growth of horizontal fractures. Physical elements such as overburden thickness, uniaxial compressive strength, rock interface position, and mining height were all shown to be related to whether horizontal cracks were present or absent. Zhao et al. [16–18] carried out a numerical analysis of the gas extraction performance of coal seams with different formations and developed a dual-system pressure decay model that can be described at different permeabilities. Tian et al., He et al., and Yang et al. [19–21] conducted indoor triaxial compression tests on rock samples under different mining disturbances, proved that mining disturbance changed the fracture height of overburden, and used numerical modeling to examine the impact of mining height and interval thickness on the height at which the fracture zone develops. Du and Gao, Ren and Wang, and Wang et al. [22–24] estimated the height and development law of the caving zone by analyzing the overburden caving characteristics of a completely automated acting surface using theoretical calculation methods, physical similarity experiments, numerical simulations, and field verification. Li et al. [25] used the DC resistivity method to dynamically monitor the fracture zone, established the numerical model and physical model of mining overburden, and assessed the height at which the fracture zone was developing before using field verification to ascertain the WCFZ’s elevation. To determine the WCFZ’s elevation in the overburdened rock of a coal mine, Zhang et al., Chai et al., and He et al. [26–28] developed a mechanical analysis model and investigated the correlations between the height of the working face, the thickness of the safety coal pillar left behind, the filling rate and other factors, and the height of WCFZ. Through progressively more in-depth theoretical derivation, physical experiments and engineering practice, many researchers have achieved a series of crucial research results on the development principles, influencing factors and height prediction of WCFZ.

The main methods to study the WCFZ are the empirical formula, the physical similarity simulation, the numerical simulation, and the field measurement [29]. Based on the key strata theory, the overburden’s key strata are harder and stronger than other layers, and when they fracture or collapse, the soft rock above them also collapses [1]. Because of this, determining the working face placement of the key strata of the overburden is the key at determining the WCFZ’s development height. First, depending on where key strata are located, this article determines the WCFZ development height in the 215102 working face of the Yue Nan coal mine. Then, verifies the precision of theoretical calculation of WCFZ development height combined with statistical formula. Finally, determines the caving zone and fracture zone’s height throughout development, respectively, through the chart of the plastic zone, the stress distribution, and the displacement variation in numerical simulation.

#### 2. Engineering Outline

The area of the Yue Nan coal mine is approximately 11.66km^{2}, the production scale is 1.2 million tons annually, and the minable coal seams in the minefield are coal seams 3, 9, and 15. From 2017 to 2020, the gas grade of the Yue Nan coal mine was identified as high gas mine. The No. 15 coal seam is now being mined at the Yue Nan coal mine. The No. 15 coal seam’s 215101 working face in mining area 1 has already been mined, and the 215102 working face is now being mined. The 215102 working face was identified as the subject of this article.

The working roadway is vertically positioned to the east of the return air roadway. The horizontal return air roadway, horizontal track roadway, and horizontal belt roadway are to the west of the workings, the north is 215101 working face, and the south is 215103 working face. The stopping line is 50 m east of the return air roadway. The general layout of the working face is rectangular in shape, with a total length of 500 m in the strike direction and 180 m in the tangential direction. The exact detailed arrangement of the working face is illustrated in Figure 1. The coal seam thickness varies from 2.50 to 4.51 m with an average thickness of 3.26 m. The structure of the coal seam is straightforward, and there are typically 0–4 layers of gangue. The entire region has a steady, mineable coal seam. The rock properties of the coal seam top are dark gray limestone, containing biological fossil debris, chert blocks, and calcite veinlets. The bottom is composed of mudstone, aluminum mudstone, and sandy mudstone, with plant fossils and star-shaped pyrite. The rock properties and thickness of the top and bottom slabs are shown in Figure 2. To control the roof, the 215102 working face is arranged in a long wall type, fully mechanized mining, one-time full height mining method, and full caving method.

#### 3. Theoretical Analysis of WCFZ Height

With the gradual progress of underground coal seam mining activities, roof behind goaf loses support of coal seam, destroying the stress balance of the overburden. The rock formation above the workings starts to move, deform, and break down due to the simultaneous action of its gravity and the overlying rock formation. The roof at the center of the goaf first reaches the tension strength limit and fails tensile, whereas the roof near the coal wall support area fails shear due to shear stress. With the increasing distance of the goaf, the rock layer above the working face moves and deform integrally from low to high. At the same time, the overburden deformation and failure display zonal variations, and the caving zone, fracture zone, and continuous deformation zone form gradually from low to high, as illustrated in Figure 3.

There are many layers of varying lithology and hardness in the rock strata above the working face, and one or more of these layers is of vital importance to the movement of the overlying rock strata as a whole or in part is known as the key strata. The former can be called the main key strata, and the latter can be called the sub key strata [8]. The destruction of the key strata will lead to the destruction of all or part of the overburden. The determination of the key strata is critical in estimating the height of the working face’s WCFZ. Xu and Qian proposed that the placement of the key strata has a substantial impact on the WCFZ’s height, carried out in-depth research on its law, and described the technique for determining the WCFZ’s development height depending on where the key strata are located [9].

##### 3.1. Determination of Overburden Hard Rock Strata

According to the theory of the composite beam, the load of the rock strata is calculated layer-by-layer from low to high. The assumption used in the study is that the weight on the rock layers is distributed uniformly. Based on the geometric characteristics of the rock stratum, such as elastic modulus, thickness, and unit weight calculate the weight of the rock strata itself and the weight generated by the overburden as a result of gravity, and comprehensively calculate the load on the rock strata. Assuming that layer 1 is a hard rock layer and its upper layer shows synchronous sinking with layer 1, while layer does not show synchronous sinking, then layer is the second hard rock layer. As layer and layer 1 show synchronous sinking, each rock layer forms a combination beam, and the load acting on the first rock layer by the combination beam theory can be derived. The load applied to it is [30]

In the equation, is the load of layer n strata acting on layer 1 strata, kPa; is the elastic modulus of layer strata, GPa; is the thickness of layer strata, m; is the gravitational acceleration, m/s^{2}; is the density of layer strata, kg/m^{3}.

Since layer and layer 1 of the rock, strata sink synchronously, and layer and layer 1 do not sink synchronously, layer itself bears part of the weight of the overburden. Therefore, the weight of strata layer and strata layer on strata layer 1 can be expressed by the following equation [31, 32]:

In the equation, is the load of layer strata acting on layer 1 strata, kPa; is the load of layer strata acting on layer 1 strata, kPa.

The geological and physical characteristics of the overburden on the 215102 working face are shown in Table 1. This data was determined by drilling samples from the workings in situ and sending them to the laboratory for testing of their mechanical parameters. Taking the mudstone layer as layer 1, calculate a load of each layer on layer 1 based on Equation (1) and Table 1. By substituting the mechanical parameters of layer 1 into Equation (1), it can be obtained that the self-load of layer 1 is kPa. By substituting the mechanical parameters of layer 2 into Equation (1), it can be obtained that the load of layer 2 acting on layer 1 is . A load of layer 3 acting on layer 1 is .

And so on, calculate the load acting on layer 1 for each layer separately, and when Equation (2) is satisfied, the layer can be identified as a hard rock layer. Based on Equation (1), calculate the load of the working face overburden layer on layer 1, and Figure 4 displays the outcomes.

Figure 4 shows a cyclical phenomenon of increasing and then decreasing rock loads. A load of layer numbers 2 ~ 4 is increasing. In layer number 5, the rock strata load suddenly decreases, with the size of 10.18 kPa, indicating that layer 5 is a hard rock stratum relative to layer 1. And so on, it is known that the hard rock layers are layer 5, layer 8, and layer 15 of the overburden.

In the overburden of the 215102 working face, there are three strata of hard rock, including 11.09 m siltstone, 20.05 m sandy mudstone, and 31.83 m sandy mudstone.

##### 3.2. Calculation of Overburden Breaking Span

The breaking span of the overburden is computed using material mechanics theory and the fixed beam model. When the maximum tensile stress of the hard rock layer m in the overburden layer exceeds the rock layer’s tensile strength, the rock layer breaks. Based on the theory of material mechanics, the stressed rock strata are simplified and the mechanical theory model of fixed support beam is adopted. The normal stress at any position in the structure of the basic top beam model is [33]

In the equation, is the positive stress at this position, MPa; is the bending moment of the section at this position, ; is the distance from this position to the neutral axis of the section, m; is the thickness of basic roof strata, m.

According to the fixed supported beam model’s features, the bending moment at both ends of the beam model is the greatest, and its maximum bending moment and maximum tensile stress are as follows:

The rock stratum will fracture and collapse, according to the model structure, when the tensile stress approaches the limit of the overburden. Then, the breaking span of layer m is can be expressed by the following equation:

In the equation, is the breaking span of layer m, m; is the thickness of layer m, m; is the tensile strength of layer m, MPa; is the load borne by layer m, MPa.

The breaking span of overburden is calculated based on Figure 4 and Equation (5). The calculation outcomes are displayed in Figure 5.

It can be seen from Figure 5 that the breaking span of the rock strata varies with the thickness, tensile strength, and load of the rock stratum, and the breaking span of the hard rock strata is larger than that of the general rock strata. The breaking span of siltstone with a thickness of 11.09 m is 7.77 m, that of sandy mudstone with a thickness of 20.05 m is 15.47 m, and that of sandy mudstone with a thickness of 31.83 m is 34.18 m.

##### 3.3. Determination of Overburden Key Strata and Calculation of WCFZ Height

After calculating the breaking span for each hard rock layer, the sizes were compared separately to identify the main and sub key stratum [9].

If the hard rock layer m is a key strata, it must have a breaking span that is smaller than all the hard rock layers above it, and it must satisfy the equation shown in the following [34]:

In the equation, is the breaking span of hard rock layer , ; is the breaking span of hard rock layer , .

It is assumed that layer and layer are hard layers through calculation, in which the layer is located above layer and , then both layer and layer are key strata. If , then layer is not key strata. Add the load carried by hard rock layer to hard rock layer , recalculate the breaking span of layer , and compare it with that of hard rock layer . If the recalculated breaking span of layer is less than that of layer , take .

Judge whether is true layer by layer from the lowest hard rock layer, and when , recalculate the breaking span of the layer .

Three hard rock strata layers have been found to exist in the overburden, and the breaking span from bottom to top is m, m, and m, judge the three key strata according to the process in Figure 6.

Because ; therefore, the three hard rock layers are broken progressively from bottom to top, and they all overburden key strata. Among them, layer 5 siltstone and layer 8 sandy mudstone are sub key strata, and layer 15 sandy mudstone is the main key strata.

Relevant experimental research indicates that the position of the key stratum has a considerable impact on the WCFZ’s height. If the distance between it and the height of the coal seam does not exceed the critical value of (7 ~ 10) M (M is the mining height of the coal seam), the key strata will break and form a fracture, and the upper rock strata will break and pass through to form a WCFZ [7]. If the spacing exceeds this critical value, then the height of the WCFZ is equal to the distance from the base of the key stratum to the working face that is greater than and most closely approximates the critical value. The overall thickness of the overburden is equal to the height of the WCFZ if the separation between the main key strata and the working face does not exceed the critical value of (7 ~ 10) M.

The mining height m is put into the calculation for the critical height by the coal seam mining data. It is evident that the critical height range is 28~40 m, which is far greater than the range from the sub key strata siltstone layer closest to the working face by 15.16 m and far less than the range from the main key strata to the working face by 87.35 m. Therefore, the key strata and their upper overburden are broken and connected to form a WCFZ, and the WCFZ’s highest growth height is 87.35 m.

##### 3.4. Empirical Formulae Validate Theoretical Calculations

To explore the correctness of the theoretical calculation results, the empirical formula is employed to determine the three overburden zones. The empirical equation may be used to determine the height of its WCFZ, and it is known that the 215102 working face belongs to the hard coal seam.

The calculation equation of caving zone is [6]

In the equation, is the height of caving zone, m; is the coal seam mining height, m.

The calculation equation of fracture zone is [6]

In the equation, is the height of fracture zone, m; is the coal seam mining height, m.

Substituting the mining height of coal seam m into Equation (7), the height range of the caving zone is 13.89~18.89 m. Substituting Equations (8) and (9), the height of the fracture zone is 49.92~67.72 m and 70 m, respectively.

The empirical formula’s calculation outcomes correspond well with the theoretical analysis results as mentioned, proving the precision of the analysis findings.

#### 4. Numerical Simulation of WCFZ Height

In order to study the development height of the water-conducting fracture zone at the 215102 working face of Yue Nan coal mine and its evolution law, numerical simulations were carried out using the 3-Dimension Distinct Element Code (3DEC) software, which is mainly applicable to the study of the mechanism of deformation and damage phenomena due to discontinuous interfaces.

As illustrated in Figure 7, the strata model of the coal seam is built in keeping with the actual geological conditions of the Yue Nan. The model size m (), the working face size is 300 m × 180 m × 4 m (), and located in the center of the model. The surroundings of the model are confined during numerical simulation, and a load with a 360-meter buried depth is imposed at the top of the model. At the bottom, the Z-directional displacement is also constrained. The calculation model adopts a linear elastic model, and the joint surface contact adopts the Coulomb sliding model. The detailed properties of each rock formation are shown in Table 1.

##### 4.1. Analysis of Overburden Plastic Zone Simulation Results

When using 3DEC software for simulation calculation, with the working face mining, the initial rock stress balance is damaged. The overburden will sink and crack under the force of gravity, and the ceiling of the goaf will likewise break and fall to regain equilibrium. The cloud map of the plastic zone can directly reflect the damage to the overburden after coal mining, and the overburden’s crack and collapse are what results in fracture development. Therefore, the distribution characteristics of the plastic zone and the development law of WCFZ correspond to each other. As a result, by scrutinizing the particular range of plastic zone, it is possible to ascertain the distribution features and development height of WCFZ.

Figure 8 depicts the features of the plastic zone’s partition at working face distances of 80 m, 120 m, 160 m, and 200 m. There are differences in the stress states reflected by various colors. Shear corresponds to shear stress, tension corresponds to tensile stress, is now to reflect the current damage in this area, is previous to reflect the previous damage in this area, and the blue area corresponds to no damage.

**(a) Push to 80 m**

**(b) Push to 120 m**

**(c) Push to 160 m**

**(d) Push to 200 m**

As seen in Figure 8, since the stope stress will be redistributed after mining the coal seam, the damage focusing on shear stress appears in the overburden, and the damage range continues to increase due to overburden fracture and collapse. Only a restricted range of coal seam roofs and floors experience shear and tensile failure in the plastic zone when the working face is pushed to 80 meters, with shear failure predominating. The failure range and degree of the plastic zone continue to expand in line with the working face’s constant development. The mining condition is sufficient when the working face is pushed 200 meters, but under the influence of the main key state, the height of the plastic zone of the overburden remains constant in comparison to when it is pushed 160 meters, indicating that the damage range of the plastic zone has expanded to its maximum. At this time, the failure area’s maximum height in the plastic zone is approximately 15.96 meters, and the failure area’s maximum height is 82.39 meters, indicating that the caving zone’s height is 15.96 meters and the fracture zone’s height is 82.39 meters.

##### 4.2. Analysis of Overburden Stress Simulation Results

The overburden is continually subjected to changing disturbance stress as the working face advances, and the rock strata fracture when the stress exceeds its elastic-plastic bearing capacity limit, resulting in cracks. Studying the dynamic process of overburden stress change during coal seam mining can therefore indirectly reflect the status of the overburden at the time, and then analyze how the overburden stability is affected by mining. As the working face continues to advance, the self-stress balance of overburden stope stress gradually changes to the state of regional pressure concentration and release. Figure 9 depicts the stope stress as the working face advances to 80 m, 120 m, 160 m, and 200 m. The stress variation graphs for pushing the face to 80 m, 120 m, 160 m and 200 m are shown in Figure 10.

**(a) Push to 80 m**

**(b) Push to 120 m**

**(c) Push to 160 m**

**(d) Push to 200 m**

**(a) Push to 80 m**

**(b) Push to 120 m**

**(c) Push to 160 m**

**(d) Push to 200 m**

As can be seen in Figures 9 and 10, with the gradual advance of the working face, a pressure relief area symmetrically distributed with the central axis of the goaf appears directly above the goaf and the floor, which expresses that the rock stratum in this area has collapsed or separated fractures. The pressure relief area shows the characteristics of high around and low in the middle. The features of high stress surrounding and low stress in the center progressively vanish when the separation between the roof and the working face increases, and the stress contour eventually becomes oval. When the working face advanced 160 m, the overburden had already started to collapse and a stress peak had appeared in the middle of the overburden in the mining area, at which point the stress peak was 11.79 MPa. When it advances to 200 m, the overburden has fully collapsed, and the development of the WCFZ is stable. The peak stress in the middle of the overburden rock in the extraction area has reached its maximum value and is close to the peak stress in the rock on both sides of the working face, where the peak stress is 20.96 MPa. The height of the stress core above the goaf is approximately 13.65 meters as the working face progresses to 200 meters, and the height of the pressure relief region is approximately 77.24 meters. This demonstrates that the maximum height of the fracture zone is about 77.24 m and the maximum height of the roof caving zone is approximately 13.65 m.

##### 4.3. Analysis of Overburden Displacement Simulation Results

The stress equilibrium of the overburden in goaf is ruined by coal seam extraction. Due to the combined force of one’s weight and the upper load, overburden breaks collapses, bends, and sinks until it contracts and compacts with the floor of goaf to reach a new equilibrium state. To study the movement and deformation of overburden after mining in 215102 working face, the change of vertical displacement will be analyzed according to the overburden, which indirectly reflects the evolution process of WCFZ.

The change of stope vertical displacement when the advancing working face reaches 80 m, 120 m, 160 m, and 200 m is shown in Figure 11. Figure 12 depicts the change in overburden displacement above the goaf as the working face advances to distances of 80, 120, 160, and 200 meters.

**(a) Push to 80 m**

**(b) Push to 120 m**

**(c) Push to 160 m**

**(d) Push to 200 m**

As the working face moves forward, as seen in Figures 11 and 12, the subsidence value of overburden in goaf gradually increases. The sinking value of overburden in goaf steadily rises and reaches a maximum of roughly 3.89 m when the working face is pushed to 160 m. The sinking value of the overburden in the goaf keeps rising as the working face is pushed out to 200 m. At this time, the maximum subsidence value is about 3.91 m, which is close to the mining height of 4 m, which proves that the mining has been stable at this time. The subsidence and collapse of the overburden will compact the goaf, and the development of WCFZ is stable. At this time, the goaf can be divided into three areas, of which 0 ~ 10.08 m is situated in the area of greatest subsidence, which is the range of the caving zone; 10.08~82.69 m is located in the area with a large subsidence value, which is the range of fracture zone; 82.69 m and above are located in the area of general subsidence value, which is the scope of the continuous deformation zone. Therefore, the fracture zone has a maximum height of around 82.69 m, whereas the caving zone’s maximum height is approximately 10.08 m.

#### 5. Discussion

The ceiling collapse behind the working face causes a significant number of cracks in the WCFZ, which alters the flow state and storage condition of the initial gas, during coal mine operations. The gas will drift upward along the caving zone in the separation zone and accumulate a large amount of gas in the fracture zone. Affected by the air leakage in the goaf, the change of atmospheric pressure, or the collapse of gangue in the goaf, the gas in the goaf flows into the coal mining face or production roadway, which affects the normal production and even leads to major accidents [35]. High level boreholes are widely used in the management of gas in the mining area, and the location of the extraction layer is particularly critical. High-level boreholes are mainly based on the “three zones” distribution of the overlying strata, and should be located in the lower and middle fracture zones where gas gushes out more intensively from adjacent layers. To avoid mine gas disasters and increase mine gas utilization rate, it is crucial to investigate WCFZ’s height in the overburden.

The aforementioned study shows that the results of the theoretical calculations based on the key strata are less different from the numerical simulations, but more different from the results of the empirical formulae. This is because the empirical formula is obtained by counting the height of the WCFZ of the coal mine for many years, which is only used as a reference. The height of the overburden WCFZ shall be determined according to the geological mining conditions of the mining area and the analysis of measured data. Therefore, the theoretical calculation and numerical simulation results based on actual geological and mechanical parameters are more accurate and have more reference value.

#### 6. Conclusions

(1)There are three key strata in the overburden of 215102 working face, including 11.09 m sub key strata siltstone, 20.05 m sub key strata sandy mudstone, and 31.83 m main key strata sandy mudstone. The findings of the theoretical calculations indicate that the WCFZ’s maximum development height is 87.35 m(2)The results of the empirical formula show that the height of the caving zone in 215102 working face is 13.89~18.89 m, and the height of the fracture zone is 49.92~67.72 m and 70 m(3)According to the simulation findings of the plastic zone in the numerical simulation of working face 215102, the caving zone’s maximum height is about 15.96 m, and the fracture zone’s maximum height is approximately 82.39 m. According to the findings of the stress simulation, the caving zone’s maximum height is about 13.65 m, and the fracture zone’s maximum height is approximately 77.24 m. According to the findings of the displacement simulation, the caving zone’s maximum height is about 10.08 m, while the fracture zone’s maximum height is approximately 82.69 m(4)It is ultimately found that the height of the overburden caving zone of 215102 working face is 14.02 m and the height of the fracture zone is 76.42 m by averaging the results of many calculation techniques

#### Data Availability

The data that support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

No potential conflict of interest was reported by the authors.

#### Acknowledgments

This research was supported by the National Natural Science Foundation of China (Nos. 51874122 and 52074105), the Key R & D and Extension Projects of Henan Province (Nos. 202102310223 and 222102320017).