Study on Water Inﬂow Variation Law of No.1 Shaft Auxiliary Shaft in HighLiGongshan Based on Dual Medium Model

: When the fracture is not developed and the connectivity is poor, the original single medium simulation cannot meet the accuracy requirements. Now, the seepage simulation of fractured rock mass has gradually developed from equivalent continuous medium to dual medium and multiple medium. However, it is still difﬁcult to establish the connection between a discrete fracture network model and a continuous medium model, which makes it difﬁcult to simulate the inﬂuence of fracture location on the seepage ﬁeld of rock mass. As the excavation direction of the shaft is vertically downward, the surrounding strata are symmetrical around the plane of the shaft axis, which is different from the horizontal tunnel. Taking the auxiliary shaft of the No.1 Shaft in HighLiGongshan as the engineering background, combined with Monte Carlo methods and DFN generator built in FLAC3D5.01, a discrete fracture network is generated. Based on the dual medium theory, MIDAS is used to optimize the modeling of each fracture group. At the same time, the concept of “Fracture Weakening area” is introduced, and the simulation is carried out based on a ﬂuid–solid coupling method. It is found that the simulation effect is close to the reality. The water inﬂow increases with the increase of shaft excavation depth, and the water inﬂow at the end of excavation is nearly three times larger than the initial value. Combined with Legendre equation, a new analytical formula of water inﬂow prediction is proposed. It is found that this analytical formula is more sensitive to permeability and has a greater safety reserve.


Introduction
A proper mathematical model must be established for fracture simulation. The commonly used mathematical models include the equivalent porous media model, the discrete fracture network model and the dual medium model. Long j.c.s., Junyi Gao, Hongjiang Chen et al. have studied the seepage-stress coupling and the seepage-heat coupling of fractured rock mass by using the equivalent porous media model [1][2][3][4][5][6][7][8][9]. The results show that the equivalent porous media model is more suitable for the fractured rock mass with the representative elementary volume (REV), but not for the large fracture with low density, and it is unable to simulate the preferential flow characteristics of fracture. Jinli Wang studied the fractured rock mass by using the discrete fracture network model [10], and found that the simulation effect of this method is better when the fracture connectivity is good and the fracture data is clear, but there are difficulties in data acquisition and connectivity evaluation. The dual media model can be divided into a narrow dual media model and a broad dual media model. Dong Yang and Chen He have studied the generalized double medium model and developed the corresponding software [11,12], but the modeling method and the establishment of the equivalent formula at the intersection are difficult, as it is not widely used. Guanglin Du, Enzhi Wang, J. rutqvist., etc., adopted the narrow dual medium model [13][14][15][16][17][18][19][20][21][22][23], established the seepage continuity equation through the principle of flow equivalence between fractures and pores, and carried out the corresponding research on the seepage model of fractured rock mass, and found that the simulation is difficult, but the effect is good.
The traditional double medium network simulation mostly uses the simple superposition of two kinds of media, ignoring the influence of fracture location, which has certain errors. With the development of computer technology, people gradually begin to study the construction of three-dimensional fracture networks. Compared with a 2D network, a 3D network considers the strike and dip angle of fractures, which makes the research more practical [24]. Now, the most commonly used method of 3D fracture network construction is the Monte Carlo method. With the deepening of research, it was proved that this method can well simulate the actual fracture network [25][26][27][28][29][30][31]. In addition, there are the LHS [32] method, random polygon method [33], spatial analysis method of deteriorated rock mass [34], artificial neural network method [35,36], etc.
The HighLiGongshan fractures belong to a steep and large fracture, so it is not suitable to adopt the theory of equivalent continuous medium. In addition, it is found that the connectivity of the fracture group is not good according to the spot investigation data. There are a large number of closed fractures, so the discrete fracture network model is also not suitable. In this paper, based on the dual medium theory, combined with the spot measured fracture group data, the corresponding data are generated by Monte Carlo methods. At the same time, the original fracture data are further discretized with the help of a FLAC3D5.01 built-in DFN plug-in to obtain the actual fracture network. The double medium theory is used for numerical simulation, and a water inflow prediction formula based on the general solution of Legendre equation is derived. It is found that the prediction result is more in line with the reality.
Considering that the distribution of fracture strike conforms to Gaussian distribution, there is symmetry to a certain extent. In addition, the shaft support is symmetrical about the axis in the horizontal direction, and it is excavated symmetrically in sections in the axial direction. Different from the tunnel, the stratum around the shaft is a vertical and downward semi-wireless space body with the surface as the interface, and the stratum around the shaft is symmetrical about the axis of the shaft. The research based on the research background of No.1 Shaft auxiliary shaft has certain symmetrical reference meanings for the subsequent construction of the No.2 Shaft nearby. The research ideas of this paper are as follows in Figure 1.

Engineering Background
The auxiliary tunnel of HighLiGongshan tunnel is set as 'one through horizontal guide + one inclined shaft + two vertical shafts', in which the auxiliary shaft of one vertical shaft is 764.74 m deep and the inner diameter of the auxiliary shaft is 5.0 m. After the completion of the auxiliary shaft, it is used for feeding, fresh air, drainage, moving personnel in and out, and also as a safety exit. At about 19:00 on 15 January 2018, the depth of the auxiliary shaft was 630.3 m. In the process of mucking out and leveling, the side of the shaft on the left side of the line suddenly collapsed, forming a 'bowl shaped' gushing point. The water inflow was about 300 m 3 /h, causing well flooding. See Figure 2a,b for the spot pictures. On 19 January 2018, Yun Gui railway company held an expert meeting. Through analysis, it was concluded that the complex hydrogeological conditions around the shaft, the sudden water inflow on the side of the shaft and the inability to accurately predict the water inflow were the main reasons for this well flooding. The water inflow of on-site engineering mainly included water gushing from the shaft wall and working face, as shown in Figure 2c,d.

Engineering Background
The auxiliary tunnel of HighLiGongshan tunnel is set as 'one through horizontal guide + one inclined shaft + two vertical shafts', in which the auxiliary shaft of one vertical shaft is 764.74 m deep and the inner diameter of the auxiliary shaft is 5.0 m. After the completion of the auxiliary shaft, it is used for feeding, fresh air, drainage, moving personnel in and out, and also as a safety exit. At about 19:00 on 15 January 2018, the depth of the auxiliary shaft was 630.3 m. In the process of mucking out and leveling, the side of the shaft on the left side of the line suddenly collapsed, forming a 'bowl shaped' gushing point. The water inflow was about 300 m 3 /h, causing well flooding. See Figure 2a,b for the spot pictures. On 19 January 2018, Yun Gui railway company held an expert meeting. Through analysis, it was concluded that the complex hydrogeological conditions around the shaft, the sudden water inflow on the side of the shaft and the inability to accurately predict the water inflow were the main reasons for this well flooding. The water inflow of on-site engineering mainly included water gushing from the shaft wall and working face, as shown in Figure 2c,d.
Due to the complex geological conditions, water inrush often occurs, and Party A decided to carry out the corresponding research. This research belongs to one of the sub topics. This research has carried on the detailed simulation and the corresponding theoretical formula derivation for the shaft water inflow problem. It is found that the power of pumping equipment should be further strengthened. Figure 3 shows the proposed drainage equipment.

Geological Conditions
The stratum exposed by the borehole is granite with an exposed thickness of 771.43 m. The lithology is mainly composed of light gray green, light gray, dark gray and light gray black granite. According to the exposed stratum thickness and core, the rock cores are dense, massive, hard, complete, locally broken, closed to micro tensional fracture, medium coarse grain metamorphic structure, grain size ranging from 0.4 to 3 cm, mainly composed  Due to the complex geological conditions, water inrush often occurs, and Party A decided to carry out the corresponding research. This research belongs to one of the sub topics. This research has carried on the detailed simulation and the corresponding theoretical formula derivation for the shaft water inflow problem. It is found that the power of pumping equipment should be further strengthened. Figure 3 shows the proposed drainage equipment.

Geological Conditions
The stratum exposed by the borehole is granite with an exposed thickness of 771.43 m. The lithology is mainly composed of light gray green, light gray, dark gray and light gray black granite. According to the exposed stratum thickness and core, the rock cores are dense, massive, hard, complete, locally broken, closed to micro tensional fracture, medium coarse grain metamorphic structure, grain size ranging from 0.4 to 3 cm, mainly composed of quartz and mica, locally weakly weathered, and locally vertical fracture the  Due to the complex geological conditions, water inrush often occurs, and Party A decided to carry out the corresponding research. This research belongs to one of the sub topics. This research has carried on the detailed simulation and the corresponding theoretical formula derivation for the shaft water inflow problem. It is found that the power of pumping equipment should be further strengthened. Figure 3 shows the proposed drainage equipment.

Geological Conditions
The stratum exposed by the borehole is granite with an exposed thickness of 771.43 m. The lithology is mainly composed of light gray green, light gray, dark gray and light gray black granite. According to the exposed stratum thickness and core, the rock cores are dense, massive, hard, complete, locally broken, closed to micro tensional fracture, medium coarse grain metamorphic structure, grain size ranging from 0.4 to 3 cm, mainly composed of quartz and mica, locally weakly weathered, and locally vertical fracture the Buck et Blasthole Figure 3. Drainage equipment.

Hydrologic Condition
According to the survey data of the No.1 Shaft supplementary exploration hole of HighLiGongshan tunnel, the logging contents include conventional logging (apparent resistivity, spontaneous potential, natural gamma, gamma gamma), inclinometer and temperature logging. Four aquifers are determined from the logging curve, which are 154. 50

Establishment of Discrete Fracture Network
In rock mechanics, the Monte Carlo method is often used to generate a random fracture network with the same distribution as the actual rock fracture in the statistical sense. The main process is: (1) collecting rock fracture parameters; (2) statistical analysis of fracture data; (3) generating random fracture network. The fractures of different depth exposed by grouting excavation in HighLiGongshan are shown in Figure 4.
According to the survey data of the No.1 Shaft supplementary exploration hole of HighLiGongshan tunnel, the logging contents include conventional logging (apparent resistivity, spontaneous potential, natural gamma, gamma gamma), inclinometer and temperature logging. Four aquifers are determined from the logging curve, which are 154. 50

Establishment of Discrete Fracture Network
In rock mechanics, the Monte Carlo method is often used to generate a random fracture network with the same distribution as the actual rock fracture in the statistical sense. The main process is: (1) collecting rock fracture parameters; (2) statistical analysis of fracture data; (3) generating random fracture network. The fractures of different depth exposed by grouting excavation in HighLiGongshan are shown in Figure 4.  In view of the bad spot environment, the data of the fracture group are mainly obtained through a segmented drilling experiment, with the drilling depth of 771 m and 66 sections. According to the recent research results and engineering practice on the geometric distribution of rock fractures, the geometric parameters of rock fractures generally follow one or several probability density distributions. Generally speaking, the location of In view of the bad spot environment, the data of the fracture group are mainly obtained through a segmented drilling experiment, with the drilling depth of 771 m and 66 sections. According to the recent research results and engineering practice on the geometric distribution of rock fractures, the geometric parameters of rock fractures generally follow one or several probability density distributions. Generally speaking, the location of fracture center points obeys uniform distribution; the trace length obeys negative exponential distribution or lognormal distribution; the occurrence usually obeys univariate or bivariate Fisher distribution, bivariate normal distribution or lognormal distribution, etc. Combined with the actual engineering background of HighLiGongshan, it is considered that the engineering fracture components are divided into three groups, the center point conforms to the uniform distribution, the trace length conforms to the normal distribution, the occurrence conforms to the normal distribution, and the density of the three groups is similar. The statistical fracture data are shown in Table 1. The (0,1) random number of uniform distribution is obtained by using the linear congruence method. The geometric characteristics of fractures can be obtained through the obtained characteristic function. For convenience, 30 fractures are selected for each group. In order to ensure the randomness of fractures, the fracture template data obtained by Monte Carlo methods are input into a FLAC3D5.01 built-in DFN program for further discretization, and the fracture network generated is shown in Figure 5.
conforms to the uniform distribution, the trace length conforms to the normal distribution, the occurrence conforms to the normal distribution, and the density of the three groups is similar. The statistical fracture data are shown in Table 1. The (0,1) random number of uniform distribution is obtained by using the linear congruence method. The geometric characteristics of fractures can be obtained through the obtained characteristic function. For convenience, 30 fractures are selected for each group. In order to ensure the randomness of fractures, the fracture template data obtained by Monte Carlo methods are input into a FLAC3D5.01 built-in DFN program for further discretization, and the fracture network generated is shown in Figure 5.

Establishment of Dual Medium Model
The concept of "fracture weakening zone" is introduced. The original fracture is simulated by the "fracture weakening zone". It is considered that the part of rock mass closely connected with long and large fractures is the "weakening zone", affected by fractures. The fluid flow in this area is dominated by fracture water, and the flow rate is approximately equal to the flow rate of fracture water. This part of rock mass is used to simulate the original fractures. In order to prove the rationality of this hypothesis, two models are adopted, one is simulated by the actual fracture width, and the other is simulated by two groups of elements contacting the fracture surface, as shown in Figure 6.

Establishment of Dual Medium Model
The concept of "fracture weakening zone" is introduced. The original fracture is simulated by the "fracture weakening zone". It is considered that the part of rock mass closely connected with long and large fractures is the "weakening zone", affected by fractures. The fluid flow in this area is dominated by fracture water, and the flow rate is approximately equal to the flow rate of fracture water. This part of rock mass is used to simulate the original fractures. In order to prove the rationality of this hypothesis, two models are adopted, one is simulated by the actual fracture width, and the other is simulated by two groups of elements contacting the fracture surface, as shown in Figure 6.
Three measuring points are taken on the fracture surface to monitor the change of pore water pressure during excavation. The seepage vector diagram is shown in Figure 7a, and the monitoring data of the measuring points are shown in Figure 7b,c.
Through Figure 7a, it is found that the fluid flow in this area is dominated by fracture water, and the flow rate is approximately equal to that of fracture water. Comparing the monitoring data, it is found that there is only a certain error in the water pressure between the two methods, but the variation law is almost the same, as shown in Figure 7. Therefore, we think this method meets the accuracy requirements of the original fracture simulation.   Three measuring points are taken on the fracture surface to monitor the cha pore water pressure during excavation. The seepage vector diagram is shown in 7a, and the monitoring data of the measuring points are shown in Figure 7b,c.  Through Figure 7a, it is found that the fluid flow in this area is dominated by fracture water, and the flow rate is approximately equal to that of fracture water. Comparing the monitoring data, it is found that there is only a certain error in the water pressure between the two methods, but the variation law is almost the same, as shown in Figure 7. Therefore, we think this method meets the accuracy requirements of the original fracture simulation.
The diameter of the auxiliary shaft of the No.1 Shaft in HighLiGongshan is 5 m. Considering the effective range of construction influence, the plane size of the model is 60 m × 60 m. As the depth of the shaft is about 770 m, considering the calculation speed of the computer, only the typical joint development depth is simulated. Combined with the spot investigation data, it is concluded that the fracture characteristics of 280~400 m and 500~620 m sections are obvious, and most of them belong to the above three groups of fractures, so a large-scale finite difference model is established based on these three sections. In addition, each fracture is accurately located through the existing fracture network. As for the introduction of the fracture network, the existing discrete network data is exported, and the fracture network is simulated by using the fine modeling function of MIDAS software. Considering the computer calculation speed and the difficulty of grid division, only the representative fractures of each group are accurately simulated, and the concept of a "fracture weakening area" is used for coarse simulation of the remaining fractures. The geological parameters of each section are shown in Tables 2 and 3.  Table 3. Geological parameters of 500 m~620 m model. The diameter of the auxiliary shaft of the No.1 Shaft in HighLiGongshan is 5 m. Considering the effective range of construction influence, the plane size of the model is 60 m × 60 m. As the depth of the shaft is about 770 m, considering the calculation speed of the computer, only the typical joint development depth is simulated. Combined with the spot investigation data, it is concluded that the fracture characteristics of 280~400 m and 500~620 m sections are obvious, and most of them belong to the above three groups of fractures, so a large-scale finite difference model is established based on these three sections. In addition, each fracture is accurately located through the existing fracture network. As for the introduction of the fracture network, the existing discrete network data is exported, and the fracture network is simulated by using the fine modeling function of MIDAS software. Considering the computer calculation speed and the difficulty of grid division, only the representative fractures of each group are accurately simulated, and the concept of a "fracture weakening area" is used for coarse simulation of the remaining fractures. The geological parameters of each section are shown in Tables 2 and 3. The geological parameters of fracture group are shown in Table 4.  -10  2  40  40  40  3  --15  3  45  45  43  4  --20 Unlike rock matrix, the seepage of fracture water is anisotropic. The seepage tensor of the three groups of fractures is obtained by Formula (1).
See Table 5 for the main value of the seepage tensor. The shaft excavation method is a drilling-blasting method, the grouting radius is 0.6 R, the excavation footage is 3.6 m, the support form is anchor mesh support, and cable element simulation is used. Considering the grouting excavation and non-grouting excavation, the interaction between MIDAS and FLAC3D is used to establish a total of four models, as shown in Figure 8, which are 1.4 million units and 660,000 units, respectively. In the model, the porous medium is symmetrical.   After the in situ stress is balanced and output by Tecplot, as shown in Figure 9, it can be seen that the existence of fractures weakens the surrounding pore water pressure.  After the in situ stress is balanced and output by Tecplot, as shown in Figure 9, it can be seen that the existence of fractures weakens the surrounding pore water pressure.  After the in situ stress is balanced and output by Tecplot, as shown in Figure 9, it can be seen that the existence of fractures weakens the surrounding pore water pressure.

Construction Plan
The wireline coring drilling tool is used for the combination construction of the survey hole. The coreless drilling tool is used to drill to the complete rock, and the coring bit is used to drill to the final hole. The geophysical logging is carried out, and the pumping test is carried out one by one in the aquifer. The equipment includes a deep well submers-

Construction Plan
The wireline coring drilling tool is used for the combination construction of the survey hole. The coreless drilling tool is used to drill to the complete rock, and the coring bit is used to drill to the final hole. The geophysical logging is carried out, and the pumping test is carried out one by one in the aquifer. The equipment includes a deep well submersible pump (model: 100QJ2-150/60), which pumps water for four times. According to the pumping time, it is divided into the fourth aquifer, the third aquifer, the second aquifer and the first aquifer.

Experimental Data
Some drilling samples are shown in Figure 10. The pumping test data are shown in Table 1, which uses the maximum drawdown depth of the water level in the well only, and the hydrogeological parameters are calculated by Formula (2). The predicted value of water inflow calculated by Formula (3) is shown in Table 2.
where Q is the single well water field, K is the permeability coefficient, M is the aquifer thickness, r w is the water level in well, 0.0455 m in this calculation, S w is the average water level draw down in stable section, and R is the radius of influence.
where Q is the average water inflow of stable section, K is the permeability coefficient, M is the aquifer thickness, h is the water level in well, s w is the average water level draw down in stable section, R is the radius of influence, and H 0 is the height from static water level of the aquifer to the aquifer floor (see Tables 6 and 7).

Basic Assumption
The strata are homogeneous porous media.
(1) The seepage conforms to Darcy's law, that is, the seepage velocity in a certain direction is proportional to the hydraulic gradient.

Formula Derivation
Let there be a differential element in the seepage field, as shown in Figure 11. The side length is d x , d y , d z , and the velocity v x , v y , v z is a function of the coordinates. The flow rate into the unit is v x , v y , v z . The flow velocity of the outflow unit is:

Formula Derivation
Let there be a differential element in the seepage field, as shown in Figure 11. The side length is , , d  Figure 11. Differential element.
It is considered that the fluid in the flow channel cannot be compressed, and the inflow flow of the unit is equal to the outflow flow. Then, it can be shown that Figure 11. Differential element.
It is considered that the fluid in the flow channel cannot be compressed, and the inflow flow of the unit is equal to the outflow flow. Then, it can be shown that i.e., ∂v x ∂x + ∂v y ∂y Darcy's law is introduced: There K is the permeability coefficient and h stand for the full head Substituting Formula (5), we get The equation expressed by column coordinate: The seepage flow in the affected area of shaft excavation is considered to be symmetrical about the z axis, that is, Equation (8) is expressed as: i.e., As there is only one variable, PDE is converted to an ordinary differential equation as follows: The solution of separated variables is as follows: The integral of Equation (12) is as follows: Using spherical coordinates to express Equation (5), we can get Take h(r, θ, ϕ) = u(r)℘(θ, ϕ) = u(r)c(θ)φ(ϕ). Substituting Equation (14), after twice separating the variables, the following formula is obtained: The seepage in the area affected by shaft excavation is considered to be symmetrical around the z-axis, i.e.; ρ = m 2 = 0 , which leads to Take b = cos θ, i.e., The general solution is as follows: The general solution of h is as follows: The general series solution of Legendre equation is as follows y = (a 0 y 1 + a 1 y 2 ) (20) because of |b| = |cos θ| ≤ 1 Formula (19) has the following form: (c l r l + b l r l+1 ) · (a 0 y 1 + a 1 y 2 )(l = 0, 1, 2, 3, · · ·, n) Based on the assumption of the radius of influence of the shaft pumping, there are Outside the ball : h 1 = +∞ ∑ l = 0 (a l r l + b l r l+1 ) · P l (b)(l = 0, 1, 2, 3, · · ·, n) In the ball : h 2 = +∞ ∑ l = 0 (c l r l + d l r l+1 ) · P l (b)(l = 0, 1, 2, 3, · · ·, n). Considering the boundary conditions, the equation is established as follows (a l r l + b l r l+1 ) · P l (b)(l = 0, 1, 2, 3, · · ·, n) When r → ±∞ , then θ → 0 ⇒ b = cos θ → 1 . Therefore, replacing it with cylindrical coordinates, we can get ln r + a(l = 0, 1, 2, 3, · · ·, n).
On the ball, d l = 0 because of r = 0 . Therefore, h 2 has the following form: It is considered that if the influence radius of shaft pumping is 2s w , (2s w , H) is a common point. From this point i.e., From this formula, we can get The general solution of the original equation is Due to b = cos θ = dh dr , a Dupuit hypothesis is adopted, which leads to After separating variables, we can get The results are as follows: where Q is the well discharge per linear meter at depth H 0 , K is the permeability coefficient, s w is the average water level draw down in stable section, r 0 is the waste path of the well, and H 0 is the ground water level at rest.

The Prediction Results of This Formula Are Compared with Those of the Original Big Well Method
The results of this formula are shown in Table 8. The comparison of the calculation results of the formula in this paper and that of the big well method is shown in Figure 12.  (32) where Q is the well discharge per linear meter at depth 0 H , K is the permeability coefficient, w s is the average water level draw down in stable section, 0 r is the waste path of the well, and 0 H is the ground water level at rest.

The Prediction Results of This Formula Are Compared with Those of the Original Big Well Method
The results of this formula are shown in Table 8. The comparison of the calculation results of the formula in this paper and that of the big well method is shown in Figure 12. Through comparison, it is found that the change of permeability coefficient in the formula derived in this paper is more sensitive. In the area with small permeability coefficient, the calculation results in this paper are very close to those calculated by the formula of pressure to no pressure. In the aquifer with large permeability coefficient, the calculation result of this formula is obviously larger than that of big well method, which has a greater safety reserve. Through comparison, it is found that the change of permeability coefficient in the formula derived in this paper is more sensitive. In the area with small permeability coefficient, the calculation results in this paper are very close to those calculated by the formula of pressure to no pressure. In the aquifer with large permeability coefficient, the calculation result of this formula is obviously larger than that of big well method, which has a greater safety reserve.

Variation Law of Water Inflow in Working Face during Shaft Excavation
It can be seen from the variation law of the monitoring curve in the Figure 13 that if the joint surface is exposed during the excavation process, the water inflow will increase by multiple, and the maximum will increase to half an order of magnitude. With the advancement of the working face, after passing through the fracture area, the water inflow of the working face will gradually fall back to the normal value, and it is in direct proportion to the excavation depth of the shaft. In addition, the simulation grouting effect is obvious, the water inflow of the working face is greatly reduced, and the water inflow basically presents the trend of stable increase with the shaft excavation, and the increase range is small. Compared with the growth rate of water inflow of the 280-400 m working face, the growth rate of water inflow of the 500-620 m working face increases significantly, which reflects that, with the increase of shaft excavation depth, the change rate of water inflow of working face increases gradually.

Variation Law of Water Inflow in Working Face during Shaft Excavation
It can be seen from the variation law of the monitoring curve in the figure 13 that if the joint surface is exposed during the excavation process, the water inflow will increase by multiple, and the maximum will increase to half an order of magnitude. With the advancement of the working face, after passing through the fracture area, the water inflow of the working face will gradually fall back to the normal value, and it is in direct proportion to the excavation depth of the shaft. In addition, the simulation grouting effect is obvious, the water inflow of the working face is greatly reduced, and the water inflow basically presents the trend of stable increase with the shaft excavation, and the increase range is small. Compared with the growth rate of water inflow of the 280-400 m working face, the growth rate of water inflow of the 500-620 m working face increases significantly, which reflects that, with the increase of shaft excavation depth, the change rate of water inflow of working face increases gradually.

Variation Law of Shaft Wall Water Inflow during Shaft Excavation
From the monitoring data in Figure 14, it can be found that the water inflow of the shaft wall increases with the increase of depth basically. The deeper the excavation depth is, the greater the increase is. When the fracture area is exposed, the water inflow of shaft wall increases to a certain extent, but the amplitude is relatively smaller than that of the working face. Combined with Figure 13, it can be found that the normal water inflow of the working face within 280-400 m is similar to that of the shaft wall, and the normal water inflow of the working face within 500-620 m is relatively smaller than that of the shaft wall. With the excavation depth, the proportion of the shaft wall water inflow to the total water inflow gradually increases from 50% to about 70%, and the proportion of shaft wall water inflow to the total water inflow is proportional to the excavation depth of the shaft. Combined with the daily discharge data, this result is in line with the reality.

Variation Law of Shaft Wall Water Inflow during Shaft Excavation
From the monitoring data in Figure 14, it can be found that the water inflow of the shaft wall increases with the increase of depth basically. The deeper the excavation depth is, the greater the increase is. When the fracture area is exposed, the water inflow of shaft wall increases to a certain extent, but the amplitude is relatively smaller than that of the working face. Combined with Figure 13, it can be found that the normal water inflow of the working face within 280-400 m is similar to that of the shaft wall, and the normal water inflow of the working face within 500-620 m is relatively smaller than that of the shaft wall. With the excavation depth, the proportion of the shaft wall water inflow to the total water inflow gradually increases from 50% to about 70%, and the proportion of shaft wall water inflow to the total water inflow is proportional to the excavation depth of the shaft.

Variation Law of Total Water Inflow during Shaft Excavation
From the monitoring data in Figure 15, it can be found that grouting in the fracture area can obviously avoid the risk of flooding caused by the rapid increase of water inflow due to the exposure of joint surface. When the fracture area is exposed by excavation, the total water inflow of the shaft increases about 2-3 times, and tends to be normal after passing through the fracture area, which is proportional to the excavation depth. With the increase of the excavation depth of the shaft, the variation range of the normal total water inflow inside the shaft gradually increases, and the variation range of 500-620 m is about 1.3 times of 280-400 m. Sufficient drainage measures should be prepared at the project site.

Comparative Analysis of Simulation, Theory and Measurement.
As the simulation depth is 280-400 m and 500-620 m, the normal total water inflow values with excavation depth of 283.6 m, 341.2 m, 388 m, 503.6 m and 608 m are taken. As the grouting method has been used in the actual project, the measured value is the water inflow after grouting excavation. It can be found from Figure 16 that the grouting simulation results are close to the field measured results, but far less than the simulation results without grouting, so the simulation effect is better. The calculation result of the original pressure to non-pressure formula is too small to meet the requirement of the safety reserve. The calculation formula in this paper is close to the simulation result of no grouting, and is closer to the actual excavation result of no grouting compared with the original pressure bearing to no pressure calculation result.

Variation Law of Total Water Inflow during Shaft Excavation
From the monitoring data in Figure 15, it can be found that grouting in the fracture area can obviously avoid the risk of flooding caused by the rapid increase of water inflow due to the exposure of joint surface. When the fracture area is exposed by excavation, the total water inflow of the shaft increases about 2-3 times, and tends to be normal after passing through the fracture area, which is proportional to the excavation depth. With the increase of the excavation depth of the shaft, the variation range of the normal total water inflow inside the shaft gradually increases, and the variation range of 500-620 m is about 1.3 times of 280-400 m. Sufficient drainage measures should be prepared at the project site.

Variation Law of Total Water Inflow during Shaft Excavation
From the monitoring data in Figure 15, it can be found that grouting in the fracture area can obviously avoid the risk of flooding caused by the rapid increase of water inflow due to the exposure of joint surface. When the fracture area is exposed by excavation, the total water inflow of the shaft increases about 2-3 times, and tends to be normal after passing through the fracture area, which is proportional to the excavation depth. With the increase of the excavation depth of the shaft, the variation range of the normal total water inflow inside the shaft gradually increases, and the variation range of 500-620 m is about 1.

Comparative Analysis of Simulation, Theory and Measurement.
As the simulation depth is 280-400 m and 500-620 m, the normal total water inflow values with excavation depth of 283.6 m, 341.2 m, 388 m, 503.6 m and 608 m are taken. As the grouting method has been used in the actual project, the measured value is the water inflow after grouting excavation. It can be found from Figure 16 that the grouting simulation results are close to the field measured results, but far less than the simulation results without grouting, so the simulation effect is better. The calculation result of the original pressure to non-pressure formula is too small to meet the requirement of the safety reserve. The calculation formula in this paper is close to the simulation result of no grouting, and is closer to the actual excavation result of no grouting compared with the original pressure bearing to no pressure calculation result.

Comparative Analysis of Simulation, Theory and Measurement.
As the simulation depth is 280-400 m and 500-620 m, the normal total water inflow values with excavation depth of 283.6 m, 341.2 m, 388 m, 503.6 m and 608 m are taken. As the grouting method has been used in the actual project, the measured value is the water inflow after grouting excavation. It can be found from Figure 16 that the grouting simulation results are close to the field measured results, but far less than the simulation results without grouting, so the simulation effect is better. The calculation result of the original pressure to non-pressure formula is too small to meet the requirement of the safety reserve. The calculation formula in this paper is close to the simulation result of no grouting, and is closer to the actual excavation result of no grouting compared with the original pressure bearing to no pressure calculation result.  Figure 16. Comparison chart of water inflow prediction.

Conclusions
(1) The simulation results of the proposed seepage model of fractured rock mass are in good agreement with the reality, which can better reflect the influence of the existence of fractures on the changes of pore pressure and water inflow. The existence of fractures will weaken the pore pressure field around, which is about 0.3 times of the pore pressure field of porous media. (2) Comparing the simulation results of grouting and non-grouting, it is found that the water inflow of grouting excavation is significantly reduced compared with that of non-grouting excavation, especially when passing through the area where the fracture is located. Compared with the monitoring data of the total water inflow in different stages, it is found that the total water inflow increases sharply when passing through the fracture area, and decreases after passing. With the progress of excavation, the total water inflow presents a decreasing trend, with the increase of shaft excavation depth, the change rate of normal total water inflow in shaft increases gradually, which is consistent with the monitoring data. (3) Compared with the theoretical calculation results, the grouting simulation results are more close to the measured water inflow, and the numerical simulation method is a good means for construction in complex water rich areas. Compared with the formula of original pressure to no pressure, the water inflow prediction formula proposed in this paper is more sensitive to the permeability coefficient, and the calculated results are larger, so it has greater safety reserves. (4) With the progress of excavation, the total water inflow of shaft increases, and the difference of water inflow before and after excavation is nearly three times. Meanwhile, the proportion of shaft wall water inflow to the total water inflow gradually increases from 50% to about 70%, and the proportion of shaft wall water inflow to the total water inflow is proportional to the excavation depth of the shaft. (5) The formula presented in this paper is applicable to the engineering with symmetrical structure and horizontal stratum such as shaft. Further research is needed for tunnel engineering. At present, the prediction of steady flow water inflow is gradually improving, and the numerical simulation and analytical solution can accurately predict the project water inflow. However, it is still difficult to predict the water inflow of a kind of unsteady flow, such as water inrush.

Conclusions
(1) The simulation results of the proposed seepage model of fractured rock mass are in good agreement with the reality, which can better reflect the influence of the existence of fractures on the changes of pore pressure and water inflow. The existence of fractures will weaken the pore pressure field around, which is about 0.3 times of the pore pressure field of porous media. (2) Comparing the simulation results of grouting and non-grouting, it is found that the water inflow of grouting excavation is significantly reduced compared with that of non-grouting excavation, especially when passing through the area where the fracture is located. Compared with the monitoring data of the total water inflow in different stages, it is found that the total water inflow increases sharply when passing through the fracture area, and decreases after passing. With the progress of excavation, the total water inflow presents a decreasing trend, with the increase of shaft excavation depth, the change rate of normal total water inflow in shaft increases gradually, which is consistent with the monitoring data. (3) Compared with the theoretical calculation results, the grouting simulation results are more close to the measured water inflow, and the numerical simulation method is a good means for construction in complex water rich areas. Compared with the formula of original pressure to no pressure, the water inflow prediction formula proposed in this paper is more sensitive to the permeability coefficient, and the calculated results are larger, so it has greater safety reserves. (4) With the progress of excavation, the total water inflow of shaft increases, and the difference of water inflow before and after excavation is nearly three times. Meanwhile, the proportion of shaft wall water inflow to the total water inflow gradually increases from 50% to about 70%, and the proportion of shaft wall water inflow to the total water inflow is proportional to the excavation depth of the shaft. (5) The formula presented in this paper is applicable to the engineering with symmetrical structure and horizontal stratum such as shaft. Further research is needed for tunnel engineering. At present, the prediction of steady flow water inflow is gradually improving, and the numerical simulation and analytical solution can accurately predict the project water inflow. However, it is still difficult to predict the water inflow of a kind of unsteady flow, such as water inrush. Institutional Review Board Statement: Not applicable.