# Analysis of Anisotropic Steady and Dynamic Seepage of Subsea Tunnels

Authors: Tianjiao Zhang
Journal: Mathematical Problems in Engineering
 (2023)

---

## Abstract

Although the underwater tunnel is under a rapid development nowadays, the theoretical research methods, including reflection and angle transformation, are still in their infancy. At present, almost all the studies are based on the hypothesis that soil layers are homogeneous, but most of the soil layers are anisotropic in actual soil engineering. Therefore, an analytical derivation method of flow field in heterogeneous media based on the traditional equiangular transformation method is presented. Numerical simulation shows that with the increase of the anisotropic permeability coefficient ratio, the contour line of approximate ellipse becomes flat. On the contrary, the lateral permeability is proportional to the upgrade of the anisotropic permeability coefficient ratio and the permeability appears the same trend.


---

## Body

## 1. Introduction

In the 1840s, the first modern underwater tunnel was built across the Thames River. Since 1930s, the establishment of underwater tunnels goes to a climax. In the past half century, the whole number of underwater tunnels existing in the world was about the sum of the previous several centuries. However, compared with the fast establishment of underwater tunnels, the theoretical development of underwater tunnels is only in the initial stage. How to extract scientific problems from engineering problems and analyze them qualitatively and quantitatively is a problem worth considering, and how to hold the major conflicts is the key factor in the theoretical study of underwater tunnels.The first attempt of mankind was in 1962 [1], when Harr first proposed the application of the mirror image method in tunnel analytical calculation, transferring the method of dealing with electrostatic field to geotechnical engineering. Subsequently, the mirror image method caused waves in the academic circles. The more representative ones were Fernandez and Alvarez [2], Lei [3], and Joo and Shin [4]. The latter simplified the lining to the boundary of an equal hydraulic head, and based on the classical Goodman solution, Goodman et al. [5] studied the relationship between the hydraulic head at the lining and the water inflow of the tunnel. The mirror image method was the application of calculation models of other disciplines in geotechnical engineering, and the tunnel was equivalent to point charge, so the influence of the grouting circle and lining cannot be described [6–12]. The advantage of the mirror image method is that the symmetry of upper and lower strata is skillfully used to analyze the seepage field, but it is only suitable for a uniform media.In recent years, as an effective method to deal with semi-infinite boundary problems, the conformal mapping method of complex variable functions has gradually developed. EL Tani [13] used the conformal mapping method to solve the seepage field of the tunnel inner wall under zero water pressure. Kolymbas and Wagner [14] assumed equal liquid energy at the outside boundary of the lining, and the theoretical solution of tunnel seepage problem was obtained by the complex variable function. As a result of using the conformal mapping method, Park et al. [15] mapped the semi-infinite seepage region into different forms of the seepage area and considered several boundary conditions to obtain theoretical solutions in some special conditions. The conformal mapping method is a theoretical method of a complex variable function. However, in the face of irregular boundary, the selection of a conformal transformation function will pose a problem [16]. For example, due to the heterogeneity of geotechnical engineering, multilayer soil will contain an inner boundary. Hu et al. [17] researched the shield tunnel in heterogeneous media and found that the interface of the media had much influence. Rieckh et al. [18] obtained a tunnel boundary element model in the anisotropic multilayer media. Zhang et al. [19] researched the influence of the deformation of the multilayer media on the tunnel. Khezri et al. [20] studied the lowest supporting pressure of tunnel in the multilayer media.Because of the different shapes and sedimentation directions of soil particles, the permeability of soil often presents the characteristics of anisotropy, and the anisotropy of infiltration will change the head distribution of the flow field and then affect the tunnel construction. At present, some scholars [14, 21] have analyzed the infinite tunnel flow field, and this paper studies the semi-infinite underwater tunnel closer to the actual situation.


## 2. Research Methods

### 2.1. Physical Model

According to the law of conformal transformation, the boundary value problem of semi-infinite space can be transformed into a finite field to be solved by a complex function. First, the physical model of the tunnel is given in Figure1.Figure 1
Physical model of a subsea tunnel.Figure1 shows the entire tunnel physical model that is established to solve the problem. In this paper, for the convenience of deduction, the origin of the Cartesian coordinate system is set on the surface of the top soil. The horizontal permeability coefficient and vertical permeability coefficient of soil are kx1 and kz1, respectively. Water depth is H and tunnel buried depth is hb. The radius of the tunnel can be expressed as r0.


### 2.2. Model Assumptions

In order to simplify the problem properly, the model studied in this paper is based on the two assumptions as follows:(i)
The tunnel is circular and located in the semi-infinite space.(ii)
The seepage field is stable and conforms to Darcy’s law.


### 2.3. Formula Derivation

The permeability coefficient of soil presents a second-order tensor distribution along all directions of space, and the tensor is expressed as follows:(1)kx100kz1,wherekx1 is the horizontal permeability coefficient and kz1 is the vertical permeability coefficient. Generalized Darcy’s law can be expressed as follows:(2)∂ψx,z∂z−∂ψx,z∂z=kx100kz1⋅∂z+γwh−z/ρg+v2/2g∂x∂z+γwh−z/ρg+v2/2g∂z.Laplace equation of seepage in the anisotropic space can be expressed as follows:(3)kx1∂2ψ∂x2+kz1∂2ψ∂z2=0.In order to make the Laplace equation of seepage flow become a standard form without coefficients, the coordinate transformation is carried out according to the following rules:(4)kx1kz1=n2=δ,where δ is defined as the anisotropic permeability coefficient ratio. The coordinate transformation rules are as follows:(5)x′=x,z′=1nz.Laplace equation after the coordinate transformation is expressed as follows:(6)∂2ψ∂x′2+∂2ψ∂z′2=0.The model after the coordinate transformation is shown in Figure2.Figure 2
Calculation model after a coordinate transformation.The conformal transformation on the model that is located in theZ plane is carried out as shown in Figure 2. The conformal transformation function is as follows:(7)Z=wζ=−iA1+ζ1−ζ.Here,(8)A=kxkzhb1−α21+α2,α=kx/kzhb−kx2/kz2hb2−r02r0.Therefore, the complex plane zoneζ is shown in Figure 3.Figure 3
Calculation model after a conformal transformation.The Laplace equation of seepage can be expressed as follows:(9)∂2ψ∂ξ2+∂2ψ∂η2=0.In the planeζ, the general solution satisfies Fourier form [14]:(10)ψ=C1+C2lnϕ+∑u=1∞C3ϕu+C4ϕ−ucosuβ.We explore the correspondence between the real number plane and complex number plane:(11)x+iz=−iA1+ζ1−ζ.The following equations can be obtained:(12)x+iz1−ξ−iη=−iA1+ξ+iη,x−ξx−iηx+iz−izξ+zη=−iA−iAξ+Aη,−iA+Aξ+Aη=−iηx−z+zξ+x−ξx+zη.By comparing the real part with the imaginary part, the equations can be obtained:(13)x−ξx+zη=Aη,A+Aξ=ηx−z+zξ,ξ=x2/A−z−z−A−x2/A−z−z−A,η=x−x2/A−z−z−A/−x2/A−z−z−A⋅xA−z.We take the topmost point of the tunnel in Figure20,−kx/kzhb−r0 for research.The coordinate value of this point in the complex plane can be obtained as follows:(14)η=x−x2/A−z−z−A/−x2/A−z−z−A⋅xA−z,ξ=kx/kzhb−r0−A−kx/kzhb−r0−A,η=0.We take the rightmost point of the tunnel for research. The coordinates of the point are as follows:(15)r0,−kxkzhb.It can be obtained that in the corresponding complex plane, the values of abscissa and ordinate are as follows:(16)ξ=0,η=r0A+kx/kzhb−r02/A+kx/kzhb+kx/kzhb−A−r02/A+kx/kzhb+kx/kzhb−A⋅r0A+kx/kzhb.The elliptic equation in the complex plane can be written as follows:(17)ξ2p2+η2q2=1.The value of the major axis is as follows:(18)p=kx/kzhb−r0−A−kx/kzhb−r0−A.The value of the minor axis is as follows:(19)q=r0A+kx/kzhb−r02/A+kx/kzhb+kx/kzhb−A−r02/A+kx/kzhb+kx/kzhb−A⋅r0A+kx/kzhb.The elliptic equation of a complex plane is expressed in the polar coordinate system:(20)ϕ=q2e2cos2β+sin2β.The boundary condition of soil surface in the complex plane can be expressed as follows:(21)ψϕ=1=H+hb.Into the general solution of the Laplace equation, you can get as follows:(22)ψϕ=1=C1+∑u=1∞C3+C4cosuβ=H+hb.So,(23)C1=H+hb,C3=−C4.The water head around the tunnel ish0, so(24)ψϕ=q2e2cos2β+sin2β=h0.Bring it into the general solution of the Laplace equation of seepage flow, and the following expression can be obtained:(25)ψϕ=q2e2cos2β+sin2β=C1+C2lnq2e2cos2β+sin2β+∑u=1∞C3q2e2cos2β+sin2βu−C3q2e2cos2β+sin2β−ucosuβ=h0.So,(26)C2=h0−H−hblnq2/e2cos2β+sin2β,C3=0.Therefore, the head of anisotropic seepage in the semi-infinite space is expressed as follows:(27)ψ=H+hb+h0−H−hblnq2/e2cos2β+sin2βlnϕ.


## 2.1. Physical Model

According to the law of conformal transformation, the boundary value problem of semi-infinite space can be transformed into a finite field to be solved by a complex function. First, the physical model of the tunnel is given in Figure1.Figure 1
Physical model of a subsea tunnel.Figure1 shows the entire tunnel physical model that is established to solve the problem. In this paper, for the convenience of deduction, the origin of the Cartesian coordinate system is set on the surface of the top soil. The horizontal permeability coefficient and vertical permeability coefficient of soil are kx1 and kz1, respectively. Water depth is H and tunnel buried depth is hb. The radius of the tunnel can be expressed as r0.


## 2.2. Model Assumptions

In order to simplify the problem properly, the model studied in this paper is based on the two assumptions as follows:(i)
The tunnel is circular and located in the semi-infinite space.(ii)
The seepage field is stable and conforms to Darcy’s law.


## 2.3. Formula Derivation

The permeability coefficient of soil presents a second-order tensor distribution along all directions of space, and the tensor is expressed as follows:(1)kx100kz1,wherekx1 is the horizontal permeability coefficient and kz1 is the vertical permeability coefficient. Generalized Darcy’s law can be expressed as follows:(2)∂ψx,z∂z−∂ψx,z∂z=kx100kz1⋅∂z+γwh−z/ρg+v2/2g∂x∂z+γwh−z/ρg+v2/2g∂z.Laplace equation of seepage in the anisotropic space can be expressed as follows:(3)kx1∂2ψ∂x2+kz1∂2ψ∂z2=0.In order to make the Laplace equation of seepage flow become a standard form without coefficients, the coordinate transformation is carried out according to the following rules:(4)kx1kz1=n2=δ,where δ is defined as the anisotropic permeability coefficient ratio. The coordinate transformation rules are as follows:(5)x′=x,z′=1nz.Laplace equation after the coordinate transformation is expressed as follows:(6)∂2ψ∂x′2+∂2ψ∂z′2=0.The model after the coordinate transformation is shown in Figure2.Figure 2
Calculation model after a coordinate transformation.The conformal transformation on the model that is located in theZ plane is carried out as shown in Figure 2. The conformal transformation function is as follows:(7)Z=wζ=−iA1+ζ1−ζ.Here,(8)A=kxkzhb1−α21+α2,α=kx/kzhb−kx2/kz2hb2−r02r0.Therefore, the complex plane zoneζ is shown in Figure 3.Figure 3
Calculation model after a conformal transformation.The Laplace equation of seepage can be expressed as follows:(9)∂2ψ∂ξ2+∂2ψ∂η2=0.In the planeζ, the general solution satisfies Fourier form [14]:(10)ψ=C1+C2lnϕ+∑u=1∞C3ϕu+C4ϕ−ucosuβ.We explore the correspondence between the real number plane and complex number plane:(11)x+iz=−iA1+ζ1−ζ.The following equations can be obtained:(12)x+iz1−ξ−iη=−iA1+ξ+iη,x−ξx−iηx+iz−izξ+zη=−iA−iAξ+Aη,−iA+Aξ+Aη=−iηx−z+zξ+x−ξx+zη.By comparing the real part with the imaginary part, the equations can be obtained:(13)x−ξx+zη=Aη,A+Aξ=ηx−z+zξ,ξ=x2/A−z−z−A−x2/A−z−z−A,η=x−x2/A−z−z−A/−x2/A−z−z−A⋅xA−z.We take the topmost point of the tunnel in Figure20,−kx/kzhb−r0 for research.The coordinate value of this point in the complex plane can be obtained as follows:(14)η=x−x2/A−z−z−A/−x2/A−z−z−A⋅xA−z,ξ=kx/kzhb−r0−A−kx/kzhb−r0−A,η=0.We take the rightmost point of the tunnel for research. The coordinates of the point are as follows:(15)r0,−kxkzhb.It can be obtained that in the corresponding complex plane, the values of abscissa and ordinate are as follows:(16)ξ=0,η=r0A+kx/kzhb−r02/A+kx/kzhb+kx/kzhb−A−r02/A+kx/kzhb+kx/kzhb−A⋅r0A+kx/kzhb.The elliptic equation in the complex plane can be written as follows:(17)ξ2p2+η2q2=1.The value of the major axis is as follows:(18)p=kx/kzhb−r0−A−kx/kzhb−r0−A.The value of the minor axis is as follows:(19)q=r0A+kx/kzhb−r02/A+kx/kzhb+kx/kzhb−A−r02/A+kx/kzhb+kx/kzhb−A⋅r0A+kx/kzhb.The elliptic equation of a complex plane is expressed in the polar coordinate system:(20)ϕ=q2e2cos2β+sin2β.The boundary condition of soil surface in the complex plane can be expressed as follows:(21)ψϕ=1=H+hb.Into the general solution of the Laplace equation, you can get as follows:(22)ψϕ=1=C1+∑u=1∞C3+C4cosuβ=H+hb.So,(23)C1=H+hb,C3=−C4.The water head around the tunnel ish0, so(24)ψϕ=q2e2cos2β+sin2β=h0.Bring it into the general solution of the Laplace equation of seepage flow, and the following expression can be obtained:(25)ψϕ=q2e2cos2β+sin2β=C1+C2lnq2e2cos2β+sin2β+∑u=1∞C3q2e2cos2β+sin2βu−C3q2e2cos2β+sin2β−ucosuβ=h0.So,(26)C2=h0−H−hblnq2/e2cos2β+sin2β,C3=0.Therefore, the head of anisotropic seepage in the semi-infinite space is expressed as follows:(27)ψ=H+hb+h0−H−hblnq2/e2cos2β+sin2βlnϕ.


## 3. Validation of Numerical Simulation

In order to verify the correctness of the analytical solution, a numerical model is established by the COMSOL software. The scale of one kilometer multiplied by one kilometer is taken to simulate the semi-infinite space. The center of the tunnel is 50 m away from the soil surface, and the radius of the tunnel is 5.5 m. The water head is applied to the soil surface through the water head boundary. The water head is 40 m, the vertical permeability coefficient of the soil layer is1.0×10−7 m/s, and the horizontal permeability coefficient is 1.0×10−7 m/s, 3.0×10−7 m/s, 5.0×10−7 m/s, and 7.0×10−7 m/s. The calculation results of head distribution are shown in Figures 4–7.Figure 4
Distribution of the hydraulic head field under the condition of isotropic seepage.Figure 5
Distribution of the hydraulic head field when the permeability coefficient ratio is 3.Figure 6
Distribution of the hydraulic head field when the permeability coefficient ratio is 5.Figure 7
Distribution of the hydraulic head field when the permeability coefficient ratio is 7.It can be seen from Figure4 that when the seepage is isotropic, the contour distribution of water head is very uniform and almost round, and the water pressure is slightly higher just below the tunnel. When the seepage becomes anisotropic, the contour line of water head becomes approximately elliptical, which is due to the slightly higher water pressure at the lower part, and with the increasing anisotropy ratio of permeability coefficient, the contour line of approximately elliptical becomes more and more flat, and the transverse seepage gradually intensifies, while the longitudinal seepage gradually weakens, which shows a trend of nonlinear change.


## 4. Analysis of Velocity

Figures8–11 show the calculation results of the water flow velocity distribution under various working conditions.Figure 8
Distribution of velocity field under the condition of isotropic seepage.Figure 9
Distribution of velocity field when the permeability coefficient ratio is 3.Figure 10
Distribution of velocity field when the permeability coefficient ratio is 5.Figure 11
Distribution of velocity field when the permeability coefficient ratio is 7.When the seepage is isotropic, the velocity distribution is very uniform, showing an approximate circular distribution, while when the seepage is anisotropic, the velocity is very uneven, and the velocity on the left and right sides of the lower part of the tunnel increases sharply. In practical engineering, the number of drainage facilities can be increased to ensure the safety of the tunnel. When the anisotropy ratio of permeability coefficient increases from 3 to 5, the change is more significant. The location where the velocity increases rapidly develops towards the lower part of the tunnel, and when the anisotropy ratio of permeability coefficient increases from 5 to 7, this phenomenon gradually weakens. In practical engineering, attention should be paid to the permeability ratio of soil to make a correct judgment.When the anisotropy ratio of soil layer increases from 1 to 3, the average velocity around the hole decreases by 17.14%, but the velocity at the bottom of the hole increases by 12.15%. When the anisotropy ratio increases from 3 to 5, the average velocity around the hole increases by 5.27%, and the velocity at the bottom slightly increases by 1.15%. When the anisotropy ratio is very large, such as 7, the average velocity around the hole increases by 16.24%, and the velocity at the bottom increases by 0.03. It can be seen that in the actual tunnel construction, the waterproofing of lining should be increased according to the anisotropy ratio of soil layer. For example, when the anisotropy ratio is greater than 3, the waterproofing of lining bottom should be significantly increased. When the anisotropy ratio continues to increase, the bottom waterproofing can take into account the slight increase in cost, but the support around the tunnel should be provided well to prevent the lateral water pressure.Take a vertical line at a certain distance (35 m) from the tunnel and compare the values of water heads from top to bottom on this line under various working conditions of permeability coefficient anisotropy ratio. The comparison of the hydraulic head distribution between analytical calculation and numerical simulation is shown in Figure12.Figure 12
Comparison between analytical calculation and numerical calculation.It can be seen from Figure12 that the analytical results of all working conditions are in good agreement with the numerical results, and there is a small error at the depth of 50 m, that is, at the center of the tunnel. With the increase of anisotropy ratio, the peak values gradually decrease to 89.15%, 78.91%, and 68.39% of the original values. This is because the numerical simulation is a finite field after all, which will affect the boundary of the tunnel opening.Figure13 shows the distribution of pore water pressure around the tunnel. Curves are analytical solutions, and scatter points are numerical solutions. The numerical solution is in good agreement with the analytical solution on the whole, and the errors at 45° and 135° are relatively large, which is due to the calculation error of pore water pressure caused by “velocity concentration” similar to stress concentration twice in the numerical calculation process. With the increase of anisotropy ratio, the peak point remains the same, but the value decreases to the original 95.13%, 87.69%, and 78.51% in turn. It can be seen that the pore pressure around the tunnel presents a similar fluctuation distribution and is very uneven. With the increase of the anisotropy ratio of permeability coefficient, the pore pressure around the tunnel presents a gradual nonlinear decrease with a large difference in peak value, which should be paid special attention to in practical engineering practice.Figure 13
Distribution of pore water pressure around the subsea tunnel.


## 5. Wave Breaking Seepage

If the boundary condition is adjusted to sinusoidal wave load, the change of the total period pore water pressure around the tunnel with the angle is shown in Figure14. It can be seen that when the anisotropic permeability coefficient ratio gradually increases, the phase difference is very significant, and with the increase of the anisotropic ratio, the peak pore water pressure gradually increases by 6.14%, 15.26%, and 27.89%. In the actual construction process, it should be combined with the phase lag of the pore water pressure around the tunnel caused by the anisotropic ratio.Figure 14
Distribution of pore water pressure around the subsea tunnel.


## 6. Conclusions

Overall, when the seepage becomes anisotropic, the water pressure in the lower part of the tunnel is slightly higher than that of the upper part. The velocity of water on both sides of the lower part increases sharply. The anisotropy ratio of permeability coefficient is found to affect the transverse seepage positively, while the longitudinal seepage gradually weakens with the increase of anisotropy ratio.Thus, in the practical engineering, the number of drainage facilities is recommended to increase to ensure the safety of the tunnel. The lining waterproofing at the positions of tunnel (45° and 135°) needs to be settled in practical projects. In the actual tunnel construction, the thickness of lining waterproofing should be increased according to the anisotropy ratio of soil layer. For example, when the value of anisotropy ratio is greater than 3, the water flow rate of lining bottom will be significantly upgraded. When the anisotropy ratio increases continually, the bottom waterproofing needs to be enhanced even though it might result in a slight increase in cost. The support around the tunnel should be provided well to prevent the lateral water pressure.



---

*Source: 4345863-2023-01-25.xml*