Study on the risk of seepage ﬁ eld of Qiantang River underground space excavated in water-rich rheological rock area

: With the rapid development of underwater space construction, the theoretical research is at the initial stage. At present, the theoretical research methods mainly include mirror image method and conformal transformation method. At present, almost all the research objects in analytical research are homogeneous soil layers. However, in practical engineering, soil layers are often distributed in layers. In this study, the seepage ﬁ eld of layered soil layers is analyzed, and the analytical derivation method of seepage ﬁ eld in heterogeneous media based on traditional conformal transformation method is given, which is veri ﬁ ed by numerical simulation. The hydraulic head loss of the medium with better permeability in shallow layer is very small, and the stream-line is vertically distributed. In the weak interlayer, the hydraulic head loss is huge, and the streamline bends greatly, showing a strong horizontal distribution and converging towards the center. It can provide theoretical basis for the design of water and soil pressure in practical engineering.


Introduction
In 2180-2160 BC, Babylon built a 900 m pedestrian tunnel across the Euphrates River from the palace to Jupiter Temple. In 1843, the first underwater tunnel (pedestrian tunnel) in modern times was built on the Thames River. In 1865, the tunnel was merged into the East London Railway and converted into a railway tunnel, which is also the first underwater railway tunnel in the world. Since 1930s, the construction of underwater tunnels has developed rapidly. In the past 50 years, the total number of underwater tunnels built in the world was almost the sum of the previous hundred years [1][2][3][4][5]. However, compared with the rapid construction speed of underwater tunnel engineering, the theoretical cognition of underwater tunnel is only in the initial stage. How to reasonably extract scientific problems from the engineering problems of underwater tunnels and how to grasp the main contradictions are the key philosophies in the theoretical research of underwater tunnels.
There were also some other methods emerging. Bobet [6,7] studied the lining stress in surrounding rock with and without water, and the analytical method was pioneered to study different seepage quantities precedent of stress of lining structure. Polubarinova-Kochina [8] derived an approximate formula for the water inflow of the circular tunnel based on the assumption that the seepage field around the tunnel is a stable flow. Zhang and Franklin [9] obtained a jointed analytical solution, which was applicable to the prediction of ubiquitous inflows from a jointed rock mass. Su et al. [10] obtained semi-analytical solutions for calculating groundwater level and tunnel water inflow. Ying et al. [11] studied the influence of boundary conditions of tidal wave on seepage field. Zareifard [12] obtained the response of tunnel seepage field under the condition of considering cracked lining.
In the recent 20 years, as a very efficient method to deal with boundary value problems, the conformal transformation method of complex variable functions has gradually matured. EL Tani [13] used conformal transformation method to map the water line and the tunnel circumference into two concentric circles, and obtained the analytical solution of seepage field under the condition of zero water pressure around the tunnel circumference. Kolymbas [14] assumed equal hydraulic head outside the lining, the analytical solution of tunnel seepage field under this condition was obtained by conformal transformation method. By using conformal transformation, Park et al. [15] transformed the semi-infinite seepage field into different forms of seepage field, and assumed different boundary conditions, thus obtaining analytical solutions in some specific cases. Conformal transformation method is a mathematical method of complex variable function, but how to apply it to geotechnical engineering and how to effectively meet the characteristics of geotechnical engineering is worth considering. The remarkable characteristic of geotechnical media is heterogeneity. Hu et al. [16] studied the stability of large-diameter shield tunnel in heterogeneous soil layer, and found that the interface of soil layer had significant influence. Rieckh et al. [17] created a tunnel boundary element model in an anisotropic multilayer media. Zhang et al. [18] studied the influence of the stiffness of multi-layer soil on the tunnel lining performance. Khezri et al. [19] calculated the minimum supporting pressure of tunnel in multilayer soil by using upper bound solution. Huang et al. [20] studied the response of multi-layer soil to moving load. Guo et al. [21][22][23] studied several kinds of geotechnical engineering problems about underwater shield tunnels and deep excavations. Yuan et al. [24][25][26] studied the sustainability of the polymer SH reinforced recycled granite residual soil and obtained displacement field of layered soils surrounding laterally loaded pile based on transparent soil. At present, there are some theoretical studies [27][28][29][30][31][32][33][34][35], and there are also some practical studies [36][37][38][39][40][41][42][43]. However, the research on seepage of tunnels in multi-layer soil is still blank, and many engineering accidents are related to groundwater. At present, although some progress has been made in numerical simulation, the analytical study is particularly important. In this work, the analytical study of heterogeneous media is mainly presented, and the analytical solution is verified by numerical software. How to create a calculation method for underwater tunnels with multi-layer media has become an urgent problem to be solved. Based on the increasingly mature conformal transformation method, this study attempts to put forward a method to study the seepage field of underwater tunnels with multi-layer media.

Tunnel physical model
According to the principle of conformal transformation method, the tunnel problem in semi-infinite space in real number field is mapped to finite complex number field to be solved. First, the physical model of tunnel in semi-infinite space in real number field is given. The following is the first analysis. Figure 1 shows the entire tunnel physical model which is established to solve the problem. In this study, for the convenience of deduction, the origin of Cartesian coordinate system is set on the surface of the top soil. From the origin of coordinates down along the Y axis, the soil layers are from the first layer to the nth layer, the permeability coefficient of which varies from k s1 to k sn , respectively. Water depth is H , and tunnel buried depth is h b . The thickness of the soil changes from h 1 to − h n 1 . The inner radius of tunnel lining is r 0 , and the outer radius is r L , and the radius of grouting ring is r G .

Model assumptions
In order to simplify the problem properly, the model studied in this study is based on the following three assumptions: (i) Each layer of the soil is isotropic.
(ii) The tunnel is circular and located in the semi-infinite space. (iii) The seepage field is stable and conforms to Darcy's law.

Analysis and deduction
The problems studied in this work have various interfaces and boundaries. The interfaces between the first layer of soil and the (n − 1) layer of soil include water-soil interface and soil-soil interface, and the interfaces between the (n -1) layer of soil include soil-soil interface and soil-grouting circle interface. Therefore, it is necessary to divide the domain. The following is an improvement of the traditional conformal transformation method [13] based on seepage theory.
3.1 Deduction of seepage field of the n-layer soil (soil-soil boundary and soil-grouting circle boundary) Let the nth layer medium be a region I n , and seepage satisfies Laplace equation so the seepage control equation in the region I n is as follows: where the physical meaning of φ sn is total hydraulic head.
Total hydraulic head can be divided into pressure hydraulic head and position hydraulic head. The pressure hydraulic head in the nth layer can be expressed as the ratio of water pressure P n to water bulk density γ w . Set the center of the tunnel as zero potential energy surface, and position hydraulic head can be expressed as follows: Initial boundary condition in the region Ι n is as follows: where H sn is the hydraulic head at the top surface of the nth soil layer. If the polar coordinate system is established in the center of the tunnel, we can get where φ L is the hydraulic head of lining area and H 0 is the hydraulic head of tunnel inner wall.
where H L is the hydraulic head of lining outer wall.
where φ G is the hydraulic head of lining area. H G is the hydraulic head of the interface between grouting ring and the nth medium (the nth soil layer).
Conformal transformation of complex variable function to physical model, and the figure after mapping is shown in Figure 2.
Conformal transformation function of complex variable function is as follows: in which In the coordinate system after conformal transformation ( − ξ η coordinate system), Laplace equation can be changed as follows: In the plane ζ , the general solution satisfies Fourier form [14].
where C 1 , C 2 , C 3 , and C 4 are constants to be determined by boundary conditions. Considering the graphic changes caused by conformal mapping, we take any point on the top surface of the Corresponding to the conformal transformation function, the following equation can be obtained: By comparing the real part with the imaginary part, the following equations can be obtained: The solution of the equations is Boundary conditions can be rewritten as follows: The flow velocity is continuous at the boundary between soil layer and grouting ring.
The following equation can be obtained: The inner boundary of grouting circle satisfies the following equation: Flow rate converged to the top surface of the (n -1) layer soil is given by Flow rate converged to the top surface of the nth layer soil is given by Continuous flow at boundary is given by So the following equation can be obtained.
3.2 Deduction of seepage field of the (n -1) layer soil (soil-soil boundary) Let the (n−1) layer medium be a region − Ι n 1 , the seepage control equation in the region − Ι n 1 is as follows: where the physical meaning of ( − ) φ n s 1 is total hydraulic head.
Total hydraulic head can be divided into pressure hydraulic head and position hydraulic head. Set the center of the tunnel as zero potential energy surface, and position hydraulic head can be expressed as follows: Boundary condition: In the coordinate system after conformal transformation ( − ξ η coordinate system), Laplace equation can be changed as follows: In the plane ζ , the general solution satisfies Fourier form.
where C 5 , C 6 , C 7 , and C 8 are constants to be determined by boundary conditions. Considering the graphic changes caused by conformal mapping, we take any point on the top surface of the (n -1)layer soil for research. We select the point with the coordi- The complex number can be expressed as follows: Corresponding to the conformal transformation function By comparing the real part with the imaginary part, the following equations can be obtained: The solution of the equations is Boundary conditions can be rewritten as follows: The bottom surface of the (n -1)-layer soil is the top surface of the n-layer soil, so

[( ) ]
(47) In the nth layer medium, By analogy, in the first and second layers of medium (soil layer), the following equations hold:   (59) cos .
Therefore, the hydraulic head distribution on the top surface of the next layer of soil is as follows: Laplace equation of grouting circle area is as follows: Boundary conditions: In the coordinate system after conformal transformation ( − ξ η coordinate system), Laplace equation can be changed as follows: In the plane ζ , the general solution satisfies Fourier form. Boundary conditions can be rewritten as follows: Bringing boundary conditions into the form of general solution, we obtain In the same way, the flow equation can be obtained due to the velocity, and flow at the boundary are continuous.    (77) The following boundary conditions and equations can be obtained for lining areas: Considering lining and grouting ring, there are + n 4 8 equations and about + n 4 8 undetermined constants, so unknown quantities (C 1 ∼ + C n 4 8 ) can be determined one by one. Therefore, the solution of Laplace equation of each layer of medium can be obtained, and the seepage field for the whole system can be determined. In other words, as long as the value of n and the permeability coefficient of each layer of soil are given, the seepage field can be determined by using the obtained equations. In order to verify the correctness of the analytical solution, a numerical model is established to simulate the seepage field of the tunnel with an engineering example. In this study, the Bo'ao tunnel project of Qiantang River in Hangzhou, China is analyzed, which mainly simulates the typical geological conditions at the bottom of Qiantang River, analyzes the water inflow and seepage field of the tunnel, studies the seepage mechanism of multilayer strata in Qiantang River, and provides guiding significance for the future construction of multi-layer underwater tunnels.

Engineering case
The geological conditions of Bo'ao Tunnel crossing the Qiantang River is shown in Figure 3. The soil is distributed   in layers, and there are three main soil layers, which are sandy silt, muddy silty clay, and rounded gravel. The rest is a small amount of interlayer. The most common water depth of Qiantang River is 8 m. Shield passes through rounded gravel layer, which is medium to dense state, the particle size is 0.5-3.0 cm. It is mainly sub-round and sub-angular, with a content of 55-60%, in which particles larger than 2 cm account for about 20%, and the particles are mainly filled with sand and clay, with sand content of about 20-30%, and silt and clay content of about 10%.  And the physical and mechanical parameters of the soil layer are shown in Table 1.
The geological conditions of the Bo'ao Tunnel crossing the Qiantang River is shown in Figure 4. Key highlight issues for analysis: the water depth H of Qiantang River is 8 m. The first layer of soil at the bottom of the river is sandy silt, the permeability coefficient of which is × − 4 10 6 m/s. The thickness of the sandy silt is 7.5 m. The tunnel is located in rounded gravel with the permeability coefficient of ×

Numerical simulation
COMSOL is used to simulate the project, and the model establishment and mesh division are shown in Figure 5. In order to ensure the accuracy and verify the analytical solution in infinite domain, the size of the numerical model is taken as 1 km × 1 km , and the situation around the cave is shown in Figure 5. COMSOL's groundwater flow module is used for simulation. This work mainly studies the seepage theory of multi-layer medium in underwater tunnel, so the stress and strain of soil is not the main research point for the time being, and the soil is regarded as elastic material to model, in order to highlight the characteristics of seepage field. The bottom is set as impervious boundary, while the left and right sides are set as normal seepage boundary. Lining and grouting ring are set to be embedded. The numerical simulation results of hydraulic head distribution are shown in Figure 6.  Hydraulic head is the representation of fluid energy per unit weight. According to the numerical simulation results, it can be seen that the hydraulic head loss is obviously different with the change in the soil layer. When the water flows through the first sandy silt, the hydraulic head loss is very small, the perpendicularity of the streamline is very high, and the water flow is basically in a state of vertical decline. However, when the water flows through muddy silty clay, the hydraulic head loss shows great changes, and the muddy silty clay layer dissipates a lot of fluid energy. At the same time, the streamline direction has also changed considerably, and the streamline suddenly becomes horizontally inclined distribution, which tends to converge in the middle in this soil layer. When the streamline enters the round gravel layer, it presents concave curve distribution and gradually merges into the tunnel. In Section 5, the hydraulic head distributed along the depth from 0 to 80 m away from the tunnel are studied, and compared with the analytical solution, which are shown in Figure 7.  It can be seen that with the change in depth, the water head distribution changes to a convex belly shape, and the change is more significant near the tunnel, which shows that near the tunnel, the hydraulic head value decreases, but the hydraulic gradient gradually increases. Moreover, the numerical solution is in good agreement with the analytical results, which shows the correctness and rationality of the analytical solution.

Parametric analysis
The difference of permeability coefficient between layered media will have different effects on the flow field. First, this study compares the difference of flow field between singlelayer soil and multi-layer soil, and compares the influence of permeability coefficient difference between multi-layer strata on seepage field. In order to compare the difference  between seepage field in single-layer soil and seepage field with weakly permeable interlayer, we take single-layer soil with permeability coefficient of × − 3 10 8 m/s and seepage field with permeability coefficient of × − 3 10 8 m/s but interlayer with permeability coefficient of × − 2 10 8 m/s to study. Figures 8-10 show the difference in simulation results.
The figures show that single-layer medium and multilayer medium show great differences in both hydraulic head distribution and flow velocity. In terms of hydraulic head distribution, the single-layer medium satisfies Laplace equation as a whole, while the multi-layer medium only satisfies the governing equation in each domain, and contains a variety of boundaries, which show great distribution differences after passing through each boundary. The influence of permeability difference of multi-layer media on seepage field is explored below. The permeability coefficient of the middle interlayer is constant, and the hydraulic head distribution is studied by taking different permeability coefficient ratios of upper and lower soil layers. The results are shown in   It can be seen that when the permeability coefficients of upper and lower soil layers change, the distribution of water head along the depth is also very different, and it has little influence at a distance from the tunnel, but the water head distribution at a position close to the tunnel (10-20 m) shows great difference. When the permeability coefficients of upper and lower soil layers are quite different, there is little difference between the water head near the tunnel. When the strong permeable soil layer is mixed with the weak permeable layer, the change in water head distribution is very obvious, and the change in permeability coefficient between the upper and lower strong permeable layers will also lead to the obvious change in water head distribution. The water head change in the actual multilayer soil engineering can be optimized for water and soil pressure with reference to the results of this study.

Conclusion
(1) Among the multi-layer media, there are significant differences in hydraulic head loss in each layer of media. The hydraulic head loss of the medium with better permeability in shallow layer is very small, and the streamline is vertically distributed. In the weak interlayer, the hydraulic head loss is huge, and the streamline bends greatly, showing a strong horizontal distribution and converging towards the center. In the soil layer with better permeability at the bottom, the streamline is distributed in concave curve, and the hydraulic head decreases toward the center of the tunnel according to a certain gradient. (2) When there is a weak interlayer, the head loss of shallow soil layer is small, while the head loss of shallow soil layer of single-layer soil is large in the middle. This also accords with the fact that the single-layer medium satisfies Laplace equation as a whole, while the multilayer medium only satisfies the governing equation in each domain. (3) For the larger value of X component of velocity, singlelayer soil is concentrated on both sides of the tunnel, while multi-layer soil is concentrated on the lower part of the tunnel. For the larger value of Y component of velocity, there are more obvious differences between single-layer soil and multi-layer soil. The direction change in single-layer soil presents concave function distribution, while presents convex function distribution in multi-layer medium. (4) For the total velocity field, the maximum value of the single-layer soil is concentrated around the tunnel and the upper part of the tunnel, while the maximum value of the multi-layer soil is concentrated in the lower part of the tunnel. (5) Under the condition of obvious difference between upper and lower soil layers, there is little difference between the change in hydraulic head distribution and permeability coefficient ratio near the tunnel. (6) With the change in permeability coefficient ratio of upper and lower soil layers, the distribution of hydraulic head along the depth is also very different, which has little effect at a certain distance from the tunnel, but the distribution of hydraulic head is very different near the tunnel.
Funding information: There is no funding for this article.
Author contributions: The author confirms sole responsibility for the following: study conception and design, data collection, analysis and interpretation of results, and manuscript preparation.