A Comparative Study on Steady-State Water Inflow into a Circular Underwater Tunnel with an Excavation Damage Zone

Abstract

Excavation damage zones that occur around tunnels usually lead to a change in rock permeability, which has an impact on the water inflow into tunnels and even induces water inrush disasters. For a better understanding of the effect of the excavation damage zone, analytical solutions that consider the excavation damage zone are developed based on the review and modification of the solutions that consider linings and grouting circles. Then, both analytical solutions and the finite element method are applied to estimate the water inflow, and the results are in good agreement. The effect of the excavation damage zone on water inflow is analyzed based on an impact factor of the excavation damage zone and evaluated in a real engineering case, and the seepage-preventing effects of grouting are discussed. The results reveal that the water inflow increases with increasing permeability and thickness of the excavation damage zone and that there is a limit for the effects of the excavation damage zone. In addition, the effect is stronger for underwater tunnels with small water inflows and stabilizes gradually as the magnitude of water inflow increases. The increase in the impermeability and thickness of the lining and grouting circle can reduce the effect.

1. Introduction

As an important option of underwater transportation engineering, underwater tunnels have the advantages of insensitivity to weather, no obstruction of navigation, good seismic performance and little impact on ecosystems. Hence, a considerable number of underwater tunnels have been constructed in the past century [1,2,3]. The estimation of water inflow is a key and difficult issue for underwater tunnels due to infinite water supply and defects of water-proofing systems [4,5,6,7]. Many efforts have been devoted to proposing effective methods to predict water inflow into tunnels over the past few decades. Four kinds of methods, including analytical solutions [8,9,10,11,12,13,14], numerical simulations [15,16,17,18,19], empirical formulas [20,21] and physical experiments [22], have been applied to estimate water inflow into tunnels.

Analytical solutions have some advantages, namely, convenience, reliability, and low cost [23]. Therefore, many researchers have developed various analytical solutions to predict water inflow into tunnels. Based on the image method proposed by Harr [24], a series of analytical solutions were derived for water inflow into tunnels in an infinite or semi-infinite aquifer [10,25,26,27]. In addition, there are several analytical and semianalytical solutions that consider the influence of lining, grouting and drawdown based on a multidomain model or numerical simulations [6,8,28]. Conformal mapping is another widely used method to obtain analytical solutions by transforming a semi-infinite aquifer into an annular or rectangular domain [23,29,30,31]. This method can deal with different boundary conditions (zero water pressure and constant total head) along the tunnel circumference and has drawn much attention from researchers. Some novel analytical solutions that consider the effects of various factors, such as lining and grouting, have been proposed based on conformal mapping [32,33,34,35,36]. Compared with the analytical solution, the numerical method can provide the ability to handle more complex geological and hydrological models and boundary conditions [13,14]. A number of studies on water inflow into tunnels using numerical methods have been conducted, and the effects of various factors on water inflow have been investigated synthetically [5,7,14,37].

It is worth noting that the formation of an excavation damage zone (EDZ) is expected to occur around tunnels due to excavation, and the EDZ is characterized by increased fractures, which lead to a significant increase in rock permeability [12,37,38,39,40]. This argument is supported by several field measurements, as follows: Souley et al. [41] reported hydraulic experiments by Canada’s Underground Research Laboratory, and the results indicated that the permeability of the EDZ can be increased by approximately 5 orders of magnitude, compared with the original rock. The results of field measurements in the Opalinus clay of the Mont Terri Rock Laboratory showed that the average hydraulic conductivity of the surrounding rock within 10–20 cm behind the lining increased by 2–4 orders of magnitude [42]. Tsang et al. [40] reported an increase in the permeability of the EDZ of up to 6 orders of magnitude for a tunnel excavated in indurated clays.

However, few of the aforementioned studies have considered the negative impact of the EDZ on water inflow into tunnels in a semi-infinite aquifer. Hence, a comparative study on steady-state water inflow into a circular underwater tunnel with an EDZ is conducted based on both analytical solutions and numerical methods. First, the analytical solutions for water inflow into a circular underwater tunnel with an EDZ in a semi-infinite aquifer have been proposed, based on the review and modification of the analytical solutions for a circular tunnel with a lining and grouting circle. Subsequently, the finite element method (FEM) is applied to estimate the water inflow into a circular underwater tunnel with an EDZ, and the calculation results of both the analytical solutions and numerical method are compared. Finally, the effect of the EDZ on water inflow into the tunnel is analyzed based on an impact factor of the EDZ proposed in this study and evaluated in the real engineering case reported in the literature, and the seepage-preventing effects of grouting are discussed.

2. Analytical Solutions

A simplified model of a circular underwater tunnel with an EDZ in a semi-infinite aquifer with a horizontal water level is shown in Figure 1. Hassani et al. [18] summarized some simplified assumptions presented in previous studies on analytical solutions for steady-state water inflow into tunnels. These simplified assumptions are modified and applied in this study as follows:

The properties of the surrounding rock, lining and EDZ are homogeneous and isotropic. The water flow is steady-state, incompressible and governed by Darcy’s law. The tunnel cross-section is circular, and the water level is horizontal.

Since the media are isotropic and homogeneous, the seepage field can be described by the following Laplace equation according to Darcy’s law and mass conservation:

$$\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}=0$$

(1)

where $\varphi$ is the total head that can be expressed as follows:

$$\varphi =\frac{p}{{\gamma }_{\mathrm{w}}}+y+h$$

(2)

where p is the water pressure, ${\gamma }_{\mathrm{w}}$ is the unit weight of water, and y + h is the elevation head.

The boundary condition at the ground surface is expressed as

$${\varphi }_{(y=0)}=H+h$$

(3)

In the case of a circular underwater tunnel, two different boundary conditions along the inner circumference of the lining were used in previous studies [29,30,31], as follows:

$${\varphi }_{(r=r_0)}=h_{\mathrm{a}}$$

(4)

$${\varphi }_{(r=r_0)}=\frac{u_{\mathrm{a}}}{{\gamma }_{\mathrm{w}}}+y+h$$

(5)

where $h_{\mathrm{a}}$ is the constant total head on the inner circumference of the lining and $u_{\mathrm{a}}$ is the constant water pressure on the inner circumference of the lining.

Considering the similarity of the EDZ and grouting circle, those analytical solutions for a circular tunnel with a lining and grouting circle in a semi-infinite aquifer with a horizontal water level can be modified and used to investigate the effect of the EDZ on water inflow into the tunnel by replacing the grouting circle with the EDZ. A brief description of the related analytical solutions can be outlined as follows:

Yang et al. [8] proposed analytical solutions for steady-state water inflow into a circular tunnel with linings and grouting circles in a semi-infinite aquifer, as follows:

$$Q=2\mathsf{\pi }\frac{H+h}{\frac{1}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h}{r_{\mathrm{G}}}}$$

(6)

where Q is the steady-state water inflow into a circular tunnel, and $k_{\mathrm{R}}$, $k_{\mathrm{L}}$, $k_{\mathrm{G}}$ are the permeability coefficients of the surrounding rock, lining, and grouting circle, respectively, and $r_{\mathrm{G}}$ is the radius of the tunnel with the lining and grouting circle.

Equation (6) is developed by considering the surrounding rock as an annular zone with an external radius of h, and the total head on the inner circumference of the lining is zero.

Du et al. [33] and Pan et al. [35] presented analytical solutions for the estimation of steady-state water inflow into a deep tunnel with a lining and grouting circle in a semi-infinite aquifer by using different conformal mapping methods. The proposed equation, assuming the total head on the inner circumference of lining h_a is constant, is as follows:

$$Q=2\mathsf{\pi }\frac{H+h-h_{\mathrm{a}}}{\frac{1}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h+\sqrt{h^2-r_{\mathrm{G}}^2}}{r_{\mathrm{G}}}}$$

(7)

Ying et al. [27] adapted the image method to derive a solution with the total head along the inner circumference of the lining of zero, as follows:

$$Q=2\mathsf{\pi }\frac{H+h}{\frac{1}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{2h}{r_{\mathrm{G}}}}$$

(8)

Zhu et al. [36] used conformal mapping to study this problem with the assumption of zero water pressure on the inner circumference of the lining, and the solution is as follows:

$$Q=2\mathsf{\pi }\frac{H+h+(A-h)C}{\frac{1}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h+\sqrt{h^2-r_{\mathrm{G}}^2}}{r_{\mathrm{G}}}}$$

(9)

where $A=\sqrt{h^2-r_{\mathrm{G}}{}^2}$ and C is a constant determined by equal seepage velocity at each intersection. The solution process and expression for C are complicated and not described in detail here. The above analytical solutions and their descriptions are listed in

Table 1. Analytical solutions for a circular tunnel with lining and grouting circle in a semi-infinite aquifer with a horizontal water table.

Reference	Equation	Description
Yang et al. [8]	$Q=2\mathsf{\pi }\frac{H+h}{\frac{1}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h}{r_{\mathrm{G}}}}$	Homogeneous and isotropic media For both shallow and deep drained tunnels ${\varphi }_{(r=r_0)}=0$
Du et al. [33]; Pan et al. [35]	$Q=2\mathsf{\pi }\frac{H+h-h_{\mathrm{a}}}{\frac{1}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h+\sqrt{h^2-r_{\mathrm{G}}^2}}{r_{\mathrm{G}}}}$	Homogeneous and isotropic media Deep tunnels ${\varphi }_{(r=r_0)}=h_{\mathrm{a}}$
Ying et al. [27]	$Q=2\mathsf{\pi }\frac{H+h}{\frac{1}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{2h}{r_{\mathrm{G}}}}$	Homogeneous and isotropic media Deep drained tunnels ${\varphi }_{(r=r_0)}=0$
Zhu et al. [36]	$Q=2\mathsf{\pi }\frac{H+h+(A-h)C}{\frac{1}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h+\sqrt{h^2-r_{\mathrm{G}}^2}}{r_{\mathrm{G}}}}$	Homogeneous and isotropic media For both shallow and deep drained tunnels ${\varphi }_{(r=r_0)}=y+h$

.

Based on the analytical solutions listed in Table 1, the analytical solutions for a circular tunnel with lining and an EDZ in a semi-infinite aquifer with a horizontal water level can be obtained by replacing the grouting circle with the EDZ, as follows:

$$Q_{\mathrm{A}1}=2\mathsf{\pi }\frac{H+h-h_{\mathrm{a}}}{\frac{1}{k_{\mathrm{E}}}\ln \frac{r_{\mathrm{E}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h}{r_{\mathrm{E}}}}$$

(10)

$$Q_{\mathrm{A}2}=2\mathsf{\pi }\frac{H+h-h_{\mathrm{a}}}{\frac{1}{k_{\mathrm{E}}}\ln \frac{r_{\mathrm{E}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h+\sqrt{h^2-r_{\mathrm{E}}^2}}{r_{\mathrm{E}}}}$$

(11)

$$Q_{\mathrm{A}3}=2\mathsf{\pi }\frac{H+h-h_{\mathrm{a}}}{\frac{1}{k_{\mathrm{E}}}\ln \frac{r_{\mathrm{E}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{2h}{r_{\mathrm{E}}}}$$

(12)

$$Q_{\mathrm{A}4}=2\pi \frac{H+h+(A-h)C}{\frac{1}{k_{\mathrm{E}}}\ln \frac{r_{\mathrm{E}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h+\sqrt{h^2-r_{\mathrm{E}}^2}}{r_{\mathrm{E}}}}$$

(13)

where $k_{\mathrm{E}}$ is the permeability coefficient of EDZ.

The constant total head on the inner circumference of lining h_a is applied for Equations (10)–(12). For Equation (13), zero water pressure on the inner circumference of the lining is used, and $A=\sqrt{h^2-r_{\mathrm{E}}{}^2}$. It is worth mentioning that a constant internal water pressure that never exceeds the initial water pressure is assumed for a tunnel with lining and a grouting circle [43]. El Tani [29] considered that the constant internal water pressure can be the production of the minimal initial water pressure by a coefficient, as follows:

$${\varphi }_{(r=r_0)}=y+h+\beta (H+h-r_0)$$

(14)

where $\beta$ is a coefficient that ranges from 0 to 1.

Therefore, the new analytical solution $Q_{\mathrm{A}4}$ that corresponds to a constant internal pore pressure can be obtained based on Equations (13) and (14), as follows:

$$Q_{\mathrm{A}4}=2\mathsf{\pi }\frac{H+h+(A-h)D-\beta (H+h-r_0)}{\frac{1}{k_{\mathrm{E}}}\ln \frac{r_{\mathrm{E}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\frac{1}{k_{\mathrm{R}}}\ln \frac{h+\sqrt{h^2-r_{\mathrm{E}}^2}}{r_{\mathrm{E}}}}$$

(15)

where D is a constant determined by equal seepage velocity at each intersection. Based on the method in reference [36], D and other constants to be determined can be acquired as follows:

$$D=\frac{1}{2}\left[E+(1-E)\frac{r_0^2}{r_{\mathrm{L}}^2}+\frac{k_{\mathrm{L}}}{k_{\mathrm{E}}}\left(E-(1-E)\frac{r_0^2}{r_{\mathrm{L}}^2}\right)\right]+\frac{1}{2}\left[E+(1-E)\frac{r_0^2}{r_{\mathrm{L}}^2}-\frac{k_{\mathrm{L}}}{k_{\mathrm{E}}}\left(E-(1-E)\frac{r_0^2}{r_{\mathrm{L}}^2}\right)\right]\frac{r_{\mathrm{L}}^2}{r_{\mathrm{E}}^2}$$

(16)

$$E=\frac{N_1+N_2{N}_3-N_4}{N_5-N_2{N}_6}$$

(17)

$$N_1=\frac{H+h-\beta (H+h-r_0)}{\frac{1}{k_{\mathrm{E}}}\ln \frac{r_{\mathrm{E}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}-\frac{1}{k_{\mathrm{R}}}\ln \alpha }(\frac{1}{r_{\mathrm{E}}}+\frac{\rho }{\alpha })$$

(18)

$$N_2=\frac{A-h}{\frac{1}{k_{\mathrm{E}}}\ln \frac{r_{\mathrm{E}}}{r_{\mathrm{L}}}+\frac{1}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}-\frac{1}{k_{\mathrm{R}}}\ln \alpha }(\frac{1}{r_{\mathrm{E}}}+\frac{\rho }{\alpha })-2k_{\mathrm{R}}A{\displaystyle \sum _{n=1}^{\infty }\frac{n{\alpha }^{n-1}({\alpha }^{2n}+1)}{{\alpha }^{2n}-1}}\rho$$

(19)

$$N_3=\frac{1}{2}\left(\frac{r_0^2}{r_{\mathrm{L}}^2}\frac{k_{\mathrm{E}}-k_{\mathrm{L}}}{k_{\mathrm{E}}}+\frac{r_0^2}{r_{\mathrm{E}}^2}\frac{k_{\mathrm{E}}+k_{\mathrm{L}}}{k_{\mathrm{E}}}\right)$$

(20)

$$N_4=\frac{1}{2}\left(\frac{r_0^2}{r_{\mathrm{L}}^2}\left(k_{\mathrm{E}}-k_{\mathrm{L}}\right)-\frac{r_0^2}{r_{\mathrm{E}}^2}\left(k_{\mathrm{E}}+k_{\mathrm{L}}\right)\right)$$

(21)

$$N_5=\frac{1}{2}\left(\frac{r_{\mathrm{L}}^2-r_0^2}{r_{\mathrm{L}}^2}{k}_{\mathrm{E}}+\frac{r_{\mathrm{L}}^2+r_0^2}{r_{\mathrm{L}}^2}{k}_{\mathrm{L}}-\frac{r_{\mathrm{L}}^2-r_0^2}{r_{\mathrm{E}}^2}{k}_{\mathrm{E}}+\frac{r_{\mathrm{L}}^2+r_0^2}{r_{\mathrm{E}}^2}{k}_{\mathrm{L}}\right)$$

(22)

$$N_6=\frac{1}{2}\left(\frac{r_{\mathrm{L}}^2-r_0^2}{r_{\mathrm{L}}^2}+\frac{r_{\mathrm{L}}^2+r_0^2}{r_{\mathrm{L}}^2}\frac{k_{\mathrm{L}}}{k_{\mathrm{E}}}+\frac{r_{\mathrm{L}}^2-r_0^2}{r_{\mathrm{E}}^2}-\frac{r_{\mathrm{L}}^2+r_0^2}{r_{\mathrm{E}}^2}\frac{k_{\mathrm{L}}}{k_{\mathrm{E}}}\right)$$

(23)

$$\alpha =\frac{r_{\mathrm{E}}}{h+A}$$

(24)

$$\rho =-\frac{2A}{{\left(h+r_{\mathrm{E}}+A\right)}^2}$$

(25)

Compared with $Q_{\mathrm{A}1}$, $Q_{\mathrm{A}2}$, and $Q_{\mathrm{A}3}$, the boundary conditions along the inner circumference of the lining of $Q_{\mathrm{A}4}$ are more realistic. To reveal the difference among these analytical solutions above, the relative error based on $Q_{\mathrm{A}4}$ is calculated by the following equation:

$${\delta }_i=\frac{Q_{\mathrm{A}i}-Q_{\mathrm{A}4}}{Q_{\mathrm{A}4}}\times 100\%$$

(26)

where i represents the number of the analytical solution and ranges from 1 to 3.

A series of computational analyses are conducted to obtain ${\delta }_i$, and the computational parameters are shown in

Table 2. Parameters of computational analyses.

Surrounding Rock		EDZ		Lining
H	0 m, 45 m	rE	4.5 m	rL	3 m
h	6~90 m	kE	1 × 10−5 m/s	r0	2.4 m
kR	1 × 10−6 m/s	-	-	kL	1 × 10−7 m/s, 1 × 10−8 m/s

.

The following assumptions are applied to make the calculation more concise:

For Equations (10)–(12), the constant total head on the inner circumference of the lining is equal to the total head at the point ($x=\pm r$, $y=-h$), namely, h_a = 0. For Equations (15), the coefficient $\beta$ is equal to zero, namely, u_a = 0.

The difference among analytical solutions with different $h/r_{\mathrm{E}}$ is shown in Figure 2. In Figure 2a,b, ${\delta }_1$ varies from 2.5% to 3%, ${\delta }_2$ varies from −0.5% to 1%, and ${\delta }_3$ varies from −1% to 0.5%, when $k_{\mathrm{R}}/k_{\mathrm{L}}=100$. In Figure 2c,d, ${\delta }_1$ varies from 13% to 24%, ${\delta }_2$ varies from 0% to 3.5%, and ${\delta }_3$ varies from −4% to 0%, when $k_{\mathrm{R}}/k_{\mathrm{L}}=10$. The above results indicate that compared with the other solutions, $Q_{\mathrm{A}1}$ overestimates the water inflow into the tunnel. ${\delta }_1$ increases as $k_{\mathrm{R}}/k_{\mathrm{L}}$ decreases and the maximum is close to 24% in this study. The underlying reason is that the surrounding rock is assumed to be an annular zone with an external radius of h to develop $Q_{\mathrm{A}1}$, which is different from the condition of a semi-infinite aquifer.

The calculation results of $Q_{\mathrm{A}2}$, $Q_{\mathrm{A}3}$ and $Q_{\mathrm{A}4}$ are quite similar, and the maximum of $\delta$ is within $\pm$4%. ${\delta }_2$ and ${\delta }_3$ are relatively large when $h/r_{\mathrm{E}}<4$, and stabilize at approximately 0% as $h/r_{\mathrm{E}}$ increases because $h+\sqrt{h^2-r_{\mathrm{E}}^2}\approx 2h$ and the water flow can be approximated as radial when $h>>r_{\mathrm{E}}$. In addition, the variation in $H/r_{\mathrm{E}}$ has little effect on $\delta$, but the variation in $k_{\mathrm{R}}/k_{\mathrm{L}}$ has a significant effect on $\delta$ according to Figure 2. When $k_{\mathrm{R}}/k_{\mathrm{L}}$ decreases from 100 to 10, the mean value of ${\delta }_1$ increases from 2.8% to 17%, and the maximum values of ${\delta }_2$ and ${\delta }_3$ increase by more than three times.

3. Numerical Method

Numerical simulation is one of the most important means to predict the water inflow into a tunnel due to its advantage in dealing with complex models and boundary conditions. A variety of numerical methods have been developed to solve problems for different media [18,19]. Considering the basic assumptions, the finite element method (FEM) is applicable to this study. Hence, the finite element commercial software package ABAQUS (ABAQUS 2021) is used to estimate the steady-state water inflow into a circular tunnel with a lining and an EDZ in a semi-infinite aquifer with a horizontal water level. According to the simplified model shown in Figure 1, a two-dimensional numerical model is created using ABAQUS/CAE 2021, as shown in Figure 3. The geometric parameters and permeabilities of the numerical model are set in the Part module and the Property module of ABAQUS/CAE 2021, according to Table 2. The other descriptions of this numerical model are as follows:

An underwater tunnel with a diameter of 6 m is excavated at h m below the ground surface from the tunnel center. The distances from the tunnel center to the bottom and both sides are (400−h) m and 500 m, respectively, which are more than 50 times the diameter. Therefore, the influence of the boundaries on water inflow can be ignored in this study. Three kinds of materials, including the surrounding rock, EDZ, and lining, are applied in this model. The Mohr-Coulomb constitutive model is adopted for the surrounding rock and EDZ, and the elastic constitutive model is adopted for the lining.

The boundary conditions of this numerical model are set in the Load module of ABAQUS/CAE 2021 and listed as follows:

The horizontal displacement of both sides is constrained, and the displacement of the bottom is constrained. The bottom boundary is impermeable, and the total heads of both sides are constant. The water pressure on the ground surface is equivalent to a uniformly distributed force (${\gamma }_{\mathrm{w}}H$). Three kinds of seepage boundary conditions along the inner circumference of lining are considered in this study, including zero total head (h_a = 0), zero water pressure ($\beta =0$), and non-zero water pressure ($\beta =0.5$).

The water inflow predicted by numerical simulations is compared with that obtained by the analytical solution, as shown in

Table 3. Water inflow by numerical method and analytical solutions.

${\mathit{k}}_{\mathbf{R}}/{\mathit{k}}_{\mathbf{L}}$	$\mathit{H}/{\mathit{r}}_{\mathbf{E}}$	$\mathit{h}/{\mathit{r}}_{\mathbf{E}}$	${\mathit{Q}}_{\mathbf{N}1}$$({\mathbf{m}}^3\cdot {\mathbf{d}}^{-1}\cdot {\mathbf{m}}^{-1})$	${\mathit{Q}}_{\mathbf{N}2(\mathit{\beta }=0)}$$({\mathbf{m}}^3\cdot {\mathbf{d}}^{-1}\cdot {\mathbf{m}}^{-1})$	${\mathit{Q}}_{\mathbf{A}4(\mathit{\beta }=0)}$$({\mathbf{m}}^3\cdot {\mathbf{d}}^{-1}\cdot {\mathbf{m}}^{-1})$	${\mathit{Q}}_{\mathbf{N}2(\mathit{\beta }=0.5)}$$({\mathbf{m}}^3\cdot {\mathbf{d}}^{-1}\cdot {\mathbf{m}}^{-1})$	${\mathit{Q}}_{\mathbf{A}4(\mathit{\beta }=0.5)}$$({\mathbf{m}}^3\cdot {\mathbf{d}}^{-1}\cdot {\mathbf{m}}^{-1})$
10	10	1.33	8.875	8.861	8.858	4.633	4.630
		2	8.156	8.149	8.112	4.252	4.233
		4	7.916	7.913	7.876	4.106	4.087
		6	8.223	8.222	8.223	4.231	4.248
		8	8.723	8.721	8.720	4.478	4.488
		10	9.294	9.293	9.277	4.770	4.762
		15	10.723	10.722	10.766	5.504	5.497
		20	12.218	12.217	12.295	6.217	6.257
	0	1.33	1.029	1.031	1.029	0.717	0.716
		2	1.377	1.352	1.346	0.854	0.850
		4	2.260	2.259	2.248	1.279	1.273
		6	3.089	3.061	3.082	1.651	1.677
		8	3.879	3.862	3.875	2.048	2.066
		10	4.641	4.646	4.638	2.446	2.442
		15	6.431	6.431	6.459	3.327	3.344
		20	8.198	8.145	8.197	4.181	4.207

. The relative error based on $Q_{\mathrm{A}4}$ is calculated by the following equation:

$${\delta }_{\mathrm{N}i}=\frac{Q_{\mathrm{N}i}-Q_{\mathrm{A}4}}{Q_{\mathrm{A}4}}\times 100\%$$

(27)

where i represents the number of the numerical simulation and ranges from 1 to 2. $Q_{\mathrm{N}1}$ represents the water inflow obtained from the numerical simulation with zero total head along the inner circumference of lining. $Q_{\mathrm{N}2}$ represents the water inflow obtained from the numerical simulation with constant water pressure along the inner circumference of the lining.

To further compare the performance of the numerical method and analytical solution, the relative error based on $Q_{\mathrm{A}4}$ for the numerical method is calculated by Equation (27), as plotted in Figure 4. In Figure 4a,b, ${\delta }_{\mathrm{N}1}$ varies from −0.8% to 2.5% and ${\delta }_{\mathrm{N}2(\beta =0)}$ varies from −0.8% to 0.6%. In Figure 4c,d, ${\delta }_{\mathrm{N}2(\beta =0.5)}$ varies from −1.6% to 0.5%. The comparison results show that $Q_{\mathrm{N}1}$ and $Q_{\mathrm{N}2(\beta =0)}$ are in very good agreement with $Q_{\mathrm{A}4(\beta =0)}$, and $Q_{\mathrm{N}2(\beta =0.5)}$ is also consistent with $Q_{\mathrm{A}4(\beta =0.5)}$, indicating that both analytical solutions and numerical methods have good performance in estimating the water inflow into a circular underwater tunnel with a lining and an EDZ in a semi-infinite aquifer.

In general, the numerical method has the advantage of dealing with the complicate geological and hydrological models with various boundary conditions, while it is more convenient and efficient to conduct parameter sensitivity analysis and summarize the law by using analytical solutions. Therefore, the effect of EDZ on water inflow into the tunnel is analyzed by the analytical solutions in the next section.

4. Discussion

4.1. Effect of EDZ on Water Inflow

As mentioned in the first section, it is very important and useful to determine the effect of the EDZ on water inflow into underwater tunnels for designing a reasonable waterproof and drainage scheme. Hence, an EDZ impact factor is proposed based on Equation (11) in this study, as follows:

$$F_{\mathrm{E}}=\frac{\ln b}{a}+\ln \left[\frac{h}{b{r}_{\mathrm{L}}}+\sqrt{{(\frac{h}{b{r}_{\mathrm{L}}})}^2-1}\right]$$

(28)

where $a=k_{\mathrm{E}}/k_{\mathrm{R}}$ and $b=r_{\mathrm{E}}/r_{\mathrm{L}}$. a is the ratio of the permeability coefficient of the EDZ to the permeability coefficient of the surrounding rock ($a\ge 1$), and b is the ratio of the radius of the EDZ to the radius of the lining ($b\ge 1$). Then, Equation (11) can be rewritten as

$$Q_{\mathrm{A}2}=\frac{2\mathsf{\pi }{k}_{\mathrm{R}}(H+h-h_0)}{F_{\mathrm{E}}+\frac{k_{\mathrm{R}}}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}}$$

(29)

For Equation (29), the ratio of water inflow into the tunnel with EDZ (Q) to the water inflow into the tunnel without EDZ ($Q_0$) can be estimated as follows:

$$\frac{Q}{Q_0}=\frac{F_{\mathrm{E}0}+\frac{k_{\mathrm{R}}}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}}{F_{\mathrm{E}}+\frac{k_{\mathrm{R}}}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}}$$

(30)

where $F_{\mathrm{E}0}=F_{\mathrm{E}(b=1)}=\ln \left[h/r_{\mathrm{L}}+\sqrt{{(h/r_{\mathrm{L}})}^2-1})\right]$.

According to Equation (30), the effect of the EDZ on water inflow into the tunnel can be evaluated based on $F_{\mathrm{E}}/F_{\mathrm{E}0}$ values. The variation in $F_{\mathrm{E}}/F_{\mathrm{E}0}$ with different parameters a and b is shown in Figure 5. Parameter a represents the ratio of permeability coefficient between the EDZ and the surrounding rock. Parameter b represents the ratio of radius between the EDZ and the lining.

As shown in Figure 5, $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases with increasing parameters a and b. $F_{\mathrm{E}}/F_{\mathrm{E}0}$ is always equal to 1, regardless of the value of a and $h/r_{\mathrm{L}}$ when $b=1$. When $b=5$ and $h/r_{\mathrm{L}}=5$, $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases from 70.21% to 0.07%, as a increases from 1 to 1000. When $b=5$ and $h/r_{\mathrm{L}}=50$, $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases from 99.95% to 65.03%, as a increases from 1 to 1000. The results above suggest that $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases with increasing a, and the variation amplitude of $F_{\mathrm{E}}/F_{\mathrm{E}0}$ increases as the EDZ becomes thicker and $h/r_{\mathrm{L}}$ decreases.

On the one hand, $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases from 100% to 70.21% as b increases from 1 to 5, when $a=1$ and $h/r_{\mathrm{L}}=5$, and it decreases from 100% to 99.95% as b increases from 1 to 5, when $a=1$ and $h/r_{\mathrm{L}}=50$. On the other hand, $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases from 100% to 0.07% as b increases from 1 to 5, when $a=1000$ and $h/r_{\mathrm{L}}=5$, and it decreases from 100% to 65.03% as b increases from 1 to 5, when $a=1000$ and $h/r_{\mathrm{L}}=50$. The results above indicate that $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases as the EDZ becomes thicker, and the variation amplitude increases as the permeability of EDZ increases and the $h/r_{\mathrm{L}}$ decreases. Therefore, the variations in EDZ parameters a and b have a greater impact on $F_{\mathrm{E}}/F_{\mathrm{E}0}$ as the tunnel is buried more shallowly.

To determine the relationship between $Q/Q_0$ and $F_{\mathrm{E}}/F_{\mathrm{E}0}$ with different $k_{\mathrm{R}}/k_{\mathrm{L}}$ and $r_{\mathrm{L}}/r_0$ values, a set of calculations are performed, and the results are shown in Figure 6 and Figure 7, respectively. Figure 6 and Figure 7 show that $Q/Q_0$ increases as $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases, namely, as a and b increase. The sensitivity of $Q/Q_0$ to $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases with increasing $k_{\mathrm{R}}/k_{\mathrm{L}}$ and $r_{\mathrm{L}}/r_0$, which indicates that the increase in the impermeability and thickness of the lining can reduce the influence of EDZ on water inflow into the tunnel. Furthermore, the increase in $Q/Q_0$ is sharp when $F_{\mathrm{E}}/F_{\mathrm{E}0}$ changes from 1 to 0.01, but the curve tends to be gentle with a further reduction in $F_{\mathrm{E}}/F_{\mathrm{E}0}$, which means that there is a limit for the effect of EDZ on water inflow into the tunnel.

According to Figure 6 and Figure 7, the $Q/Q_0$ of the deep tunnel is greater than that of shallow tunnel when $F_{\mathrm{E}}/F_{\mathrm{E}0}$ is the same. However, the opposite phenomenon is observed when the EDZ parameters are the same because $F_{\mathrm{E}}/F_{\mathrm{E}0}$ is more sensitive to the variations in EDZ parameters for shallow tunnels, as shown in

Table 4. Computational analyses for $Q/Q_0$ with different $h/r_{\mathrm{L}}$, $k_{\mathrm{R}}/k_{\mathrm{L}}$ and $r_{\mathrm{L}}/r_0$ values.

Scheme No.	$\mathit{h}/{\mathit{r}}_{\mathbf{L}}$	a	b	${\mathit{F}}_{\mathbf{E}}/{\mathit{F}}_{\mathbf{E}0}$	${\mathit{k}}_{\mathbf{R}}/{\mathit{k}}_{\mathbf{L}}$	${\mathit{r}}_{\mathbf{L}}/{\mathit{r}}_0$	$\mathit{Q}/{\mathit{Q}}_0$
1	5	1000	5	0.07%	1	1.25	11.19
	50	1000	5	65.03%	1	1.25	1.50
2	5	1000	5	0.07%	10	1.25	2.03
	50	1000	5	65.03%	10	1.25	1.31
3	5	1000	5	0.07%	1	2	4.30
	50	1000	5	65.03%	1	2	1.44

. In addition, the difference in $Q/Q_0$ between tunnels with $h/r_{\mathrm{L}}$ values of 5 and 50 decreases significantly when $k_{\mathrm{R}}/k_{\mathrm{L}}$ increases by 10 times, when comparing the data of the first group and the second group in Table 4. Similarly, the difference in $Q/Q_0$ between shallow and deep tunnels decreases with the increasing $r_{\mathrm{L}}/r_0$ values, when comparing the data of the first group and the third group in Table 4.

4.2. Effect Evaluation in Real Engineering Case

According to the above analysis, the effect of EDZ on water inflow into underwater tunnels cannot be ignored. To provide practical insight into the effect of EDZ, the water inflow into the Qingdao Jiaozhou Bay Subsea Tunnel is calculated by Equation (11), and the effect of the EDZ is evaluated by Equation (30), based on the data from the literature [8]. The parameters of the Qingdao Jiaozhou Bay Subsea Tunnel are listed in

Table 5. Parameters of Qingdao Jiaozhou Bay Subsea Tunnel.

${\mathit{r}}_0$ (m)	${\mathit{r}}_{\mathbf{L}}$ (m)	${\mathit{k}}_{\mathbf{L}}$ (m/s)	$\mathit{a}$	$\mathit{b}$
8	9	1 × 10−7	10	1.5

, and the parameters of sections and evaluation results are shown in

Table 6. Parameters of sections and evaluation results.

No.	Mileage	${\mathit{k}}_{\mathbf{R}}$ (m/s)	$\mathit{h}$ (m)	$\mathit{H}$ (m)	$\mathit{Q}$$({\mathbf{m}}^3\cdot {\mathbf{d}}^{-1}\cdot {\mathbf{m}}^{-1})$	$\mathit{Q}/{\mathit{Q}}_0$
1	ZK3 + 125–ZK3 + 325	1.74 × 10−7	27	0	1.629	1.259
2	ZK3 + 325–ZK3 + 595	8.1 × 10−7	30.3	5	6.356	1.162
3	ZK3 + 595–ZK3 + 932	4.63 × 10−7	40.6	5.3	4.906	1.162
4	ZK3 + 932–ZK4 + 102	3.47 × 10−7	46.4	6.6	4.241	1.160
5	ZK4 + 102–ZK4 + 362	4.63 × 10−7	44.9	12.4	5.860	1.154
6	ZK4 + 362–ZK4 + 422	1.5 × 10−6	44.9	16.1	13.520	1.103
7	ZK4 + 422–ZK4 + 562	4.05 × 10−7	43.7	19.8	5.915	1.161
8	ZK4 + 562–ZK4 + 627	6.94 × 10−7	43.64	21.86	9.145	1.140
9	ZK4 + 627–ZK4 + 917	3.47 × 10−7	43.24	26.26	5.742	1.166
10	ZK4 + 917–ZK4 + 977	4.63 × 10−7	43.5	29.3	7.548	1.156
11	ZK4 + 977–ZK5 + 022	4.63 × 10−7	43.86	30.64	7.697	1.156
12	ZK5 + 022–ZK5 + 212	9.26 × 10−7	44.37	33.43	13.08	1.127
13	ZK5 + 212–ZK5 + 372	3.47 × 10−7	43	37.2	6.643	1.167
14	ZK5 + 372–ZK5 + 482	1.16 × 10−7	41.7	40.8	2.631	1.193
15	ZK5 + 482–ZK6 + 082	1.16 × 10−7	40.5	42.4	2.685	1.197

. In this study, the EDZ parameters a and b are taken as 10 and 1.5, with reference to the numerical simulation in the Section 3 due to the lack of field measurement data of the EDZ.

The $Q/Q_0$ evaluated with the corresponding water inflow is plotted in Figure 8, where $Q/Q_0$ ranges from 1.103 to 1.259 in these 15 sections of the Qingdao Jiaozhou Bay Subsea Tunnel. The maximum $Q/Q_0$ occurs in section No.1, with a $Q$ value of 1.629 ${\mathrm{m}}^3\cdot {\mathrm{d}}^{-1}\cdot {\mathrm{m}}^{-1}$, and the minimum $Q/Q_0$ occurs in section No.6, with a $Q$ value of 13.52 ${\mathrm{m}}^3\cdot {\mathrm{d}}^{-1}\cdot {\mathrm{m}}^{-1}$. To further investigate the relationship, the variation in $Q/Q_0$ with the magnitude of water inflow into tunnel $Q$ is illustrated in Figure 9, which shows that $Q/Q_0$ depends on $Q$ via a power function. In Figure 9, the effect of EDZ decreases with the increase in magnitude of water inflow into the tunnel, and the rate of variation decreases gradually, indicating that the effect of EDZ is stronger for underwater tunnels with small water inflow and stabilizes as the magnitude of water inflow increases. In addition, the evaluation results, based on the Qingdao Jiaozhou Bay Subsea Tunnel, suggest that the EDZ increases the water inflow into a tunnel by at least 10%, even if the lining is applied, and it is necessary to consider the effect of EDZ carefully for estimating the water inflow into an underwater tunnel. It is worth mentioning that there is no relation between $Q/Q_0$ and the depth of water $H$ according to Equation (30), but the initial water inflow $Q_0$ increases as $H$ increases, which means that a greater water inflow $Q$ occurs in an underwater tunnel with a greater water depth.

4.3. Seepage-Preventing Effect of Grouting

Considering that the water inflow into the tunnel and the risk of water inrush increase due to the existence of the EDZ, it is necessary to take seepage-preventing measures to reduce the effect of EDZ. Grouting reinforcement during tunnel excavation can effectively improve the impermeability of the EDZ and reduce the water inflow into the tunnel. Based on the simplified model of water inflow into underwater tunnels with an EDZ (Figure 1), when the radius of the grouting circle is equal to or greater than the radius of EDZ, the analytical solution for water inflow into an underwater tunnel with a grouting circle can be obtained by replacing $k_{\mathrm{E}}$ and $r_{\mathrm{E}}$ in Equation (11) with $k_{\mathrm{G}}$ and $r_{\mathrm{G}}$. The seepage-preventing effect of grouting can be estimated by the ratio of water inflow into the tunnel with a grouting circle (Q_G) to water inflow into the tunnel with an EDZ (Q), as follows:

$$\frac{Q_{\mathrm{G}}}{Q}=\frac{\frac{k_{\mathrm{R}}}{k_{\mathrm{E}}}\ln \frac{r_{\mathrm{E}}}{r_{\mathrm{L}}}+\frac{k_{\mathrm{R}}}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\ln \frac{h+\sqrt{h^2-r_{\mathrm{E}}^2}}{r_{\mathrm{E}}}}{\frac{k_{\mathrm{R}}}{k_{\mathrm{G}}}\ln \frac{r_{\mathrm{G}}}{r_{\mathrm{L}}}+\frac{k_{\mathrm{R}}}{k_{\mathrm{L}}}\ln \frac{r_{\mathrm{L}}}{r_0}+\ln \frac{h+\sqrt{h^2-r_{\mathrm{G}}^2}}{r_{\mathrm{G}}}}$$

(31)

A series of computational analyses for $Q_{\mathrm{G}}/Q$ with different grouting parameters are performed based on Equation (31) and the input data of computational analyses are listed in

Table 7. Input data of computational analyses for $Q_{\mathrm{G}}/Q$ with different grouting parameters.

Scheme No.	${\mathit{k}}_{\mathbf{R}}/{\mathit{k}}_{\mathbf{E}}$	${\mathit{r}}_{\mathbf{E}}/{\mathit{r}}_{\mathbf{L}}$	${\mathit{k}}_{\mathbf{R}}/{\mathit{k}}_{\mathbf{G}}$	${\mathit{r}}_{\mathbf{G}}/{\mathit{r}}_{\mathbf{L}}$	${\mathit{k}}_{\mathbf{R}}/{\mathit{k}}_{\mathbf{L}}$	${\mathit{r}}_{\mathbf{L}}/{\mathit{r}}_0$	$\mathit{h}/{\mathit{r}}_{\mathbf{E}}$	$\mathit{h}/{\mathit{r}}_{\mathbf{G}}$
1	0.1	2	10	2	10	1.25	2.5	2.5
	0.1	2	100	2	10	1.25	2.5	2.5
2	0.01	2	10	2	10	1.25	2.5	2.5
	0.01	2	10	4	10	1.25	2.5	2.5
3	0.1	4	10	4	10	1.25	1.25	1.25
	0.1	2	10	4	10	1.25	1.25	1.25
	0.01	2	10	4	10	1.25	1.25	1.25

. The results of computational analyses are illustrated in Figure 10. According to the results of the first group, when $k_{\mathrm{R}}/k_{\mathrm{E}}=0.1$, $Q_{\mathrm{G}}/Q$ decreases from 67.32% to 16.95% as $k_{\mathrm{R}}/k_{\mathrm{G}}$ increases from 10 to 100, which suggests that the increase in the impermeability of the grouting circle can reduce the effect of EDZ on water inflow. Similarly, the results of the second group show that $Q_{\mathrm{G}}/Q$ decreases from 67.02% to 50.39% as $r_{\mathrm{G}}/r_{\mathrm{L}}$ increases from 2 to 4, when $k_{\mathrm{R}}/k_{\mathrm{E}}=0.01$, which means that the water inflow into the tunnel can be reduced by increasing the thickness of the grouting circle. The variation range of the permeability of the grouting circle is much larger than that of its thickness in practical engineering. Hence, reducing the water inflow by improving the impermeability of the grouting circle on the premise of $r_{\mathrm{G}}\ge r_{\mathrm{E}}$ is preferred. Moreover, the results of the third group indicate that the initial values of $k_{\mathrm{R}}/k_{\mathrm{E}}$ and $r_{\mathrm{E}}/r_{\mathrm{L}}$ have little effect on $Q_{\mathrm{G}}/Q$.

The above discussion on reducing the effect of the EDZ with grouting reinforcement is for a general model with simplified conditions. According to the previous studies [37,38,39,40,41,42], the changes in rock permeability caused by excavation damage depend on a variety of conditions, including the initial stress field, the properties of rock, the natural fracture zones and the excavation methods. Therefore, in order to put forward the specific anti-seepage measures for a real engineering case, a detailed study on how excavation damage causes the changes in rock permeability should be conducted based on field and indoor tests in the future.

5. Conclusions

In this study, both analytical solutions and numerical methods are applied to estimate the steady-state water inflow into a circular underwater tunnel with an EDZ. The comparison results show that both analytical solutions and numerical methods have good performance. An impact factor of the EDZ is proposed, and the effect of EDZ on water inflow into the tunnel is analyzed by using the analytical solution due to its convenience. Finally, a series of parameter analyses is performed, and some main conclusions are summarized as follows:

• The $F_{\mathrm{E}}/F_{\mathrm{E}0}$ ratio decreases as the permeability and thickness of the EDZ increase, and the variation amplitude increases with the decrease in $h/r_{\mathrm{L}}$. Hence, the variations in EDZ parameters a and b have a greater impact on $F_{\mathrm{E}}/F_{\mathrm{E}0}$, as the tunnel is buried more shallowly.
• The $Q/Q_0$ ratio increases as $F_{\mathrm{E}}/F_{\mathrm{E}0}$ decreases, namely, as a and b increase. However, there is a limit for the influence of the EDZ on water inflow into the tunnel, and the increase in the impermeability and thickness of the lining can reduce the effect of the EDZ on water inflow.
• The difference in $Q/Q_0$ values between deep and shallow tunnels decreases with increasing impermeability and lining thickness. Furthermore, the effect of the EDZ is stronger for underwater tunnels with small water inflows and stabilizes as the magnitude of the water inflow increases.
• Grouting is an effective measure to eliminate the effect of the EDZ on water inflow into tunnels. The increase in the impermeability and thickness of the grouting circle can reduce water inflow into the tunnel, and the initial values of $k_{\mathrm{R}}/k_{\mathrm{E}}$ and $r_{\mathrm{E}}/r_{\mathrm{L}}$ have little effect on $Q_{\mathrm{G}}/Q$. In addition, reducing the water inflow by improving the impermeability of the grouting circle on the premise of $r_{\mathrm{G}}\ge r_{\mathrm{E}}$ in practical engineering is preferred.