The Influence of Hydraulic Conditions on the Stability of Dual Circular Tunnels in Unsaturated Soils

,e performance of geotechnical structures in unsaturated soils is affected significantly by the hydraulic conditions. In the present paper, a unified computational upper bound limit analysis method is applied to study the stability of dual circular tunnels located in unsaturated soils. ,e linings are substituted by an equivalent uniform pressure exerted on the periphery of the tunnel. ,e main focus is put on the effect of ground water table and surface water infiltration on the required supporting pressure and collapse mechanisms of dual tunnels.,e critical center-to-center spacing above which the interaction of the tunnels disappears is also discussed.


Introduction
e use of underground space has becoming overwhelmingly popular with the rapid development of urbanization. Tunnels are one of the most important underground infrastructures which can significantly alleviate the traffic pressures and can also be applied in the form of pipelines and so on. e construction of dual circular tunnels is often believed to be a better option than a single large circular tunnel for practical and economic concerns. On the contrary, due to the limitation of topographical conditions, closely spaced dual tunnels have to be constructed. e interaction of the two tunnels leads to an increase of the complexity in stability analysis. In the past, a number of studies have been reported on the analysis of dual tunnels by various methods such as theoretical derivations, numerical analysis, and model tests.
e ground movements around parallel tunnels in clays were investigated by Wu and Lee [1] by a series of centrifuge model tests. e tunnels were supported by the inside compressed air pressure, and the displacement of the soil mass was continuously monitored until collapsed. Field measurements on twin tunnels excavated by the earthpressure-balanced (EPB) shields in silty soils were reported by Chen et al. [2]. Based on field data and numerical analysis, Mirhabibi and Soroush [3] discussed the effect of surface buildings on the ground settlement of twin tunnels. Reduced-scale model tests were reported by Kim et al. [4] to study the effect of tunnel proximity and alignment, liner stiffness, and soil properties on the interactions between closely spaced tunnels in clays.
Numerical methods can be employed to model complicated phenomena efficiently and flexibly. Computational limit analysis (CLA) is a direct method for prediction of the critical state of structures, which combines the flexibility of numerical methods and the rigorous property of bounding theorems. It has seen a lot of improvements in the past few decades [5][6][7][8][9]. Recently, CLA was frequently applied in the stability analysis of twin tunnels under the condition of plain strain [10][11][12][13][14][15][16], and the emphasis was mainly given to the impact of the factors such as critical tunnel spacing and novel section shapes. e stability of unsupported and uniform pressures-supported dual tunnels in soils was determined by Sahoo and Kumar [10,11] using computational upper bound limit analysis. Both purely cohesive and cohesive-frictional soils were considered and the results were normalized and presented in stability charts. Sahoo and Kumar [12] discussed the behavior of two circular parallel tunnels in clay with an overlay of sand. It is observed that the variation of the required supporting pressure in layered soils depends on the combined influence of the unit weight ratio, the soil thickness, and the friction angle of the overlying sand. Rigorous lower and upper bound solutions of the limit surcharge loading exerted on the top of dual circular tunnels in cohesive-frictional soils were obtained by Yamamoto et al. [13], which were also presented in stability charts for practical use. Rigid-block mechanisms were established and the critical spacing between the two tunnels was discussed in detail. Later, Yamamoto et al. [14] analyzed the mechanical performance of dual square tunnels in cohesive-frictional materials which subject to surcharge loading using the same numerical tools. e effect of tunnel spacing on the stability of undrained dual circular tunnels was investigated by Wilson et al. [15] using numerical limit analysis and semianalytical rigid-block techniques. Stability numbers of dual unlined elliptical tunnels in cohesive-frictional soils were obtained by Zhang et al. [16] using rigid translatory elements.
Most available studies mainly dealt with the surrounding soils as single-phased materials, and therefore, their shear strengths are kept constant. However, unsaturated soils are frequently encountered in geotechnical engineering, and a large portion of their strengths is contributed by the matric suction. e matric suction is affected by the external hydraulic conditions such as rainfall, irrigation, and ground water level fluctuation. erefore, the stability of dual tunnels in unsaturated soils will be varying with the hydraulic states of surrounding soils, which will be the main focus of the present study. Computational limit analysis is one of the most appropriate methods in geotechnical stability analysis. Sloan et al. have made important contributions to computational limit analysis in not only mathematical programming techniques but also spatial approximations and adaptive remeshing rules [5,6,8,[17][18][19][20]. Recently, Yuan and Du [21,22] incorporated the suction stress-based effective stress into the classical bounding theorems to extend the limit analysis framework to the stability analysis of unsaturated soils. One of the most important contributions in the authors' work is to consider the saturation variation of unsaturated soils by the incorporation of equivalent forces, whereas the strength parameters are specified in advance.
e proposed upper bound theory in Yuan and Du [21] will be combined with the finite element method in the present paper to make use of the seepage analysis programs of available commercial software. Parameter analysis is then conducted to discuss the effect of ground water table, the magnitude of surface water infiltration, the buried depth of the tunnels, hydraulic properties of the surrounding soils, and the tunnel spacing on the mechanical behavior of twin tunnels. e rest of this paper is organized as follows. In Section 2, the suction stress-based effective stress and the unified computational limit analysis using the finite element method and the second-order cone programming are introduced.
Hydraulic properties of the surrounding soils are given in Section 3, and the established method is verified in Section 4. Section 5 discusses in detail the effect of ground water table and surface water infiltration on the required supporting pressures, collapse mechanisms, and critical center-tocenter spacings of tunnels in three different soils. Concluding remarks are provided in Section 6.

Formulations
2.1. Basic Assumptions. In limit analysis, the deformations of the surrounding soils are assumed to be very small and all the variables are computed in their initial positions. As a result, the produced errors in the present solutions become obvious if the deformation of surrounding soils is significant. Furthermore, the present analysis aims to analyze the stability of the surrounding soils and the deformation and the collapse of the supporting structures are not considered.

Upper Bound eory for Unsaturated Soils.
CLA is a direct method for determination of limit loads based on the classical bounding theory. e tedious iterative process in finite element simulations and the experience-oriented hypotheses regarding the failure modes needed in the limit equilibrium method are no longer necessary. CLA will be an appropriate technique in the present analysis of dual tunnels in unsaturated soils if unsaturated seepage analysis is incorporated. A unified upper bound method for unsaturated soils has been recently established by Yuan and Du [21] through the introduction of the effective stress proposed by Lu et al. [23,24]: where σ ′ is the effective stress, σ is the total stress, u a is the pore air pressure, I � 1 1 0 T is an identity vector, and σ s is the suction stress. e closed form solution of the suction stress is [24] where u w is the pore water pressure and α and n are the fitting parameters of the soil water characteristic curve (SWCC) presented by Mualem [25] and van Genuchten [26]. It was shown in [21] that the shearing strength increment contributed by the matric suction can be considered by an additional work rate done by the suction stress. e soils are then dealt with as equivalent single-phased materials and the strength parameters are invariant. e following mathematical programming problem is established to obtain the limit loads of unsaturated earth structures: 2 Mathematical Problems in Engineering maximize λ, where λ is the load multiplier, _ ε is the strain rate, _ u is the velocity, g and T are the unified body force and the unified surface force specified beforehand, _ ε v is the volumetric strain rate, and f is the failure criterion expressed in terms of the effective stress. If the soils are saturated, equation (3) becomes maximize λ, It is observed that equation (4) turns into the formulation for limit analysis of saturated soils when pore water pressure is considered. erefore, unsaturated soils could be consistently analyzed by the present model.

Numerical Discretization.
e finite element method is utilized in the present paper to discretize the solution domain. A simplex strain element will be used to improve numerical efficiency and to overcome volumetric locking [7]. If all the sides are straight, the strain rate and the stress within a triangle belong to the simplex defined by their values at the corner nodes. e velocities in an element are expressed in terms of the element velocity vector _ d where _ u and _ v are the components of the velocity and N u is a coefficient matrix consisting of shape functions of the velocities. Correspondingly, the strain rate in some element is given as where L i is the area coordinate with respect to the ith vertex, N ε is a coefficient matrix consisting of shape functions of the strains, and _ ε e is the element strain rate vector. e effective stress within an element is where N σ is the coefficient matrix for the effective stresses, σ ′ e is the element effective stress vector, and σ i ′ is the effective stress at the vertexes. e strain rate within an element can be obtained by the element velocity vector as where B is the strain matrix and L is a differential operator relating the strain rate and the velocities. Introduction of where N e is the number of elements and N m is the number of element boundaries where surface tractions are specified. Equation (9) can be simplified as maximize λ where σ ′ g is the global effective stress vector, q is the global external force vector, and q s is the equivalent load vector of the suction stress. It should be mentioned that, although the stresses are employed as optimization variables, equation (10) is in fact the dual formulation of the upper bound theory. e stress-based upper bound optimization problem has been proven by Krabbenhoft et al. [27] and Makrodimopoulos and Martin [7] to be able to yield a rigorous upper bound solution.

Second-Order Cone Programming.
e soils are supposed to satisfy the Mohr-Coulomb yield criterion: where τ f and σ n ′ are the tangential and normal effective stresses on the failure surface and c ′ and φ ′ are the effective cohesion and angle of internal friction, respectively. e minus sign in front of the normal stress appears because the stress is supposed to be negative in compression. e second-order cone programming is an efficient solution technique for large-scale nonlinear optimization problems which will be utilized herein for equation (10). A modified effective stress vector s ′ is firstly introduced where σ m ′ is the mean effective stress, σ xx ′ and σ yy ′ are the components of the effective stress, s ij is the deviatoric stress, and δ ij is Kronecker's delta. For the condition of plain strain, the failure criterion is rewritten as Mathematical Problems in Engineering An auxiliary variable z is introduced to reformulate equation (13) into where a � sin φ ′ and b � c ′ cos φ ′ . e discretized unified upper bound formulation will then be summarized as maximize λ where s ′ g is the global modified effective stress vector.
Equation (15) will be solved by the MATLAB optimization tool of Mosek, which is frequently used by the researchers and has been proved efficient for large-scale problems.

Hydraulic Properties
Hydraulic properties of unsaturated soils including the permeability function and the soil water characteristic curve have to be specified firstly to determine the distribution of pore water pressure before limit analysis can be performed. e pore air pressure is assumed to be zero in all the following discussions for simplicity. Gardner's model [28] is used to define the permeability of unsaturated soils: where k w and k s are the permeabilities of unsaturated and saturated soils. Genuchten's [26] soil water characteristic curve is employed to describe the relation between the effective saturation S e and the matric suction: where α and n have the same meanings as those in equation (2).

Verification
e established numerical formulation is verified at first through the stability analysis of the unsaturated foundation shown in Figure 1. A uniform pressure is exerted on the top of the foundation. To compare the results with those in the literature [21], the ground water table is assumed to be located at 0 m, 5 m, 10 m, and 15 m below the top ground, and the analysis is performed under hydrostatic conditions. e solution domain is chosen as that in [21]. Hydraulic parameters of the foundation are taken as n � 3, α � 0.01 kPa −1 , and the gravity of the foundation soil is 20 kN/m 3 . Half of the domain is considered for the symmetry of the problem and the used mesh is also given in Figure 1, where red lines denote the fixed boundary and the green line indicates the symmetry centerline. e obtained normalized bearing capacities are compared with those in [21] in Figure 2. e close agreement between those results proves the accuracy of the present formulation, which can be therefore utilized to study the effect of hydraulic states on the stability of unsaturated soils.

Problem Definition.
e configuration of the dual tunnels studied in the present paper is given in Figure 3. e geometrical and material properties and the external loads are assumed to be invariant along the longitudinal direction of the infinite long tunnel, and therefore, the problem is solved under the condition of plain strain. e buried depth and the diameter of the tunnels are H and D, and the centerto-center distance between the two tunnels is S. No surcharge loading is applied on the top ground and a uniform supporting pressure exerted on the periphery of the tunnel is used to substitute the action of support systems. e supporting pressures of the two tunnels are supposed to be identical. Under the gravity of surrounding soils, there exists a critical pressure, below which the tunnel can no longer maintain its stability. is pressure is sometimes referred to as the loosening pressure of the tunnel and usually deemed as one of the most important parameters in tunnel designs. Obviously, CLA is very appropriate to compute the required supporting pressure and if unsaturated seepage analysis is combined, the impact of hydraulic conditions such as ground water table fluctuation and surface water infiltration can be studied efficiently. e distance between the ground water table and the top ground is H w . e water table will be placed at some typical locations, indicated as P A − P J , as shown in Figure 3

Single-Phased Dual Tunnels.
In the first step, the stability numbers of dual unlined tunnels without the incorporation of pore pressures are computed to verify the established model. e dimensionless stability number is defined as N � c lim D/c, where c lim is the maximum unit weight of the surrounding soils the tunnels can resist without support. e obtained results for the angles of internal friction of 10°, 20°, and 30°are listed in Figure 4 and compared with those in the literature. It is shown that all results agree very well and therefore the established model can be reliable for stability analysis of dual tunnels in soils. Table. During the construction and operation of tunnels, the location of the ground water table is changing all the time, leading to a variation of the hydraulic states of surrounding soils and correspondingly the variation of the loosening pressures. erefore, the range of the ground water table must be determined in design of dual tunnels. Figure 5 gives the suction stress profiles of different soils under hydrostatic conditions, where y is the distance above the ground water table. Hydraulic properties of those soils are taken as those mean values summarized by Likos et al. [29] and listed in Table 1. It is shown that the magnitude of suction stress in silt, clay, and sandy clay increases with the increase of the distance above the ground water table. e largest suction stress is seen in clay, whereas the suction stress in sand is the smallest. e differences of the suction stress of those soils become more obvious with the vertical distance. It can be seen from the right part of Figure 5 that a peak value exists in the suction stress of sand, which decreases to zero eventually. e maximum suction stress of sand is only about 500 Pa, and therefore, the effect of the hydraulic states on the stability of the tunnels is caused mainly by the variation of the unit weight of surrounding soils.

Effect of Ground Water
In the following discussions, hydraulic parameters of the surrounding soils are taken as those of silt. e effective cohesion is 10 kPa and the effective angle of internal friction is chosen to be 10°, 20°, and 30°, respectively. e porosity of those soils is 0.5, and the specific weight of soil particles is 2.7. It should be mentioned that the hydraulic properties of silt are selected only as representative values and the different sets of strength properties are used to provide guidance to practical engineering. Variations of the limit supporting pressures with tunnel spacing for H/D � 1 are given in Figure 6. A similar trend is observed of the limit supporting pressures for all ground water tables, which decrease gradually with the increase of the tunnel spacing.
ere exists a critical spacing, above which the stability of the two tunnels will be independent of each other. As the angle of internal friction is low (10°or 20°), the drop of the ground water table will result in a reduction of the critical spacing. For example, the critical spacing changes from 5D to 3.5D when the water table falls from the top ground to 40 m below it. After the friction angle attains 30°, however, an increase of the critical spacing is found with the decrease of the ground water table. e above phenomenon can be attributed to the fact that the shearing strength of unsaturated soils contributed by the matric suction is related to the angle of internal friction. Smaller supporting pressures are required to maintain the stability of the tunnels as the angle of internal friction increases, and the decreasing rate leaps with the drawdown of the ground water table.
e limit supporting pressure varies between 180 kPa and 210 kPa if the surrounding soils are all saturated, and it decreases from 100 kPa for the friction angle of 10°to −20 kPa for the friction angle of 30°when the tunnel spacing is 2D and the ground water table is placed at P G . e negative value of the limit supporting pressure means the tunnel can maintain its stability without any support measures. erefore, it can be concluded that the location of the ground water table has a large impact on the limit supporting pressure and the critical spacing of dual tunnels.
Variations of the limit supporting pressure with the ground water table for H/D � 1 are given in Figure 7. For all tunnel spacings and all angles of internal friction, the loosening pressure decreases with the fall of the ground water table, due to the contribution of the matric suction. A larger declining rate is seen when the ground water table is located above the bottom of the tunnel because the collapsed domain is confined to that region. It is interesting to see from Figure 7 that some curves coincide within a given range of ground water table and then separate when the water table changes, which clearly shows the effect of ground water table on the interaction of the dual tunnels. e obtained plastic dissipations corresponding to two cases of S/D � 2 and S/ D � 3 with the friction angle of 20°are given in Figures 8 and  9 , respectively. It can be seen from Figure 8 that although the  Mathematical Problems in Engineering location of the ground water table has a large impact on the limit supporting pressure, the collapse mechanism, however, is unrelated to the water table. e situation is different in the case of S/D � 3, where it can be noticed that, as H w /D is less than 2, the collapse domain extends to the center line of the dual tunnels, indicating their mutual influence, whereas the collapsed region is limited to the top of the tunnel when H w /D is larger than 2. erefore, the effect of the ground water table on the collapse mechanism depends on the material properties of the surrounding soils and the geometric properties of the tunnels. Variations of the limit supporting pressures with tunnel spacing and the ground water table for H/D � 3 are given in Figures 10 and 11 , and the corresponding plots for H/D � 5 are given in Figures 12 and 13 , respectively. e same form of variation is observed as that in the case of H/D � 1. e critical center-to-center tunnel spacing increases with the increase of the buried depth of the tunnels. e limit supporting pressures are still varying with tunnel spacing for S/ D � 6, φ ′ � 10°, and H/D � 3, and the steady value has not reached for S/D � 8, φ ′ � 10°and H/D � 5. erefore, the buried depth has to be considered in determination of the   Zhang et al. [16] Present paper    Mathematical Problems in Engineering critical tunnel spacing. e variation of the limit supporting pressure is more significant when the ground water table is located within the vertical domain of the tunnels, which can be attributed to the fact that the stability is mostly affected by the shearing strength of the soils near the tunnels. e effect of the water content can be negligible if the water table is below the tunnel and the distance between them is large. Furthermore, the discrepancies of the limit supporting pressure induced by tunnel spacing are more prominent when the ground water table falls below the bottom of the tunnel. erefore, it can be appropriate to say that the drawdown of the ground water table aggravates the influence of tunnel spacing.

Effect of Surface Water Infiltration.
Surface water infiltration caused by rainfall or irrigation is another ingredient that will lead to the cyclic variation of hydraulic states of the surrounding soils, especially for those shallow buried tunnels. In this section, the effect of surface water infiltration on the required supporting pressures of shallow buried tunnels will be discussed. e hydraulic properties of three soils are chosen, namely, sand, clay, and silt, as given in Table 1. Strength properties of those soils are supposed to be identical to make sure the differences in mechanical behavior of dual tunnels subject to surface water infiltration are caused only by the variation of the hydraulic states. e buried depth of the tunnels is D and the ground water table is fixed at the bottom of the tunnel, namely, the point "P D " in Figure 3, which is assumed to be independent of the surface water infiltration. A steady infiltration is considered, with the magnitude varying between 0 and 4.71 × 10 −8 m/s. For simplicity, the analysis is performed under the condition of one-dimensional vertical flow, and the analytical solution of the suction stress is given as [30] σ s � 1 α where q is the infiltration rate and c w is the unit weight of the pore water. e distribution of the suction stress of the three soils which subject to water infiltration is shown in Figure 14. It is seen that the suction stress can be reduced significantly by water infiltration, which will be uniform in space under the condition of steady infiltration. e extent of reduction is related to the type of the soil. e suction stress in clay is eliminated almost completely if the infiltration rate is −4.71 × 10 −8 m/s, 94% of its saturated permeability coefficient. e suction stress in sand even increases with the increase of the infiltration rate, although the magnitude of the suction stress is negligible in comparison with those of other soils. e increase of the suction stress in sand is due to the fact that when the sand is very dry, the suction stress can be very small and then it is increases with the saturation. Variations of the limit supporting pressures for φ ′ � 10°a re given in Figure 15, where it can be observed that the supporting pressure decreases with the increase of the tunnel spacing before reaching the steady value. e critical centerto-center spacing is reduced by the water infiltrating into the surrounding soils. e variations of the supporting pressures with the infiltration rate for φ ′ � 10°are given in Figure 16. Obviously, the required supporting pressures are increased with the intensification of the infiltration for all soils. Although the suction stress in sand improves to some extent under infiltration, the increase of the unit weight of surrounding soils deteriorates the stability of the tunnels. e limit supporting pressure grows fast at the instant of infiltration, and then, the changing rate almost reaches a constant value. e dual tunnels in different surrounding soils are affected to different degrees by water infiltration. e maximum difference of critical spacing caused by infiltration Mathematical Problems in Engineering is D for clay but 0.5 D for silt and sand; the supporting pressure increments of those soils are 34 kPa, 23 kPa, and 1.8 kPa, respectively, due to their different hydraulic behavior under infiltration, as shown in Figure 14. e limit supporting pressures for φ ′ � 20°and φ ′ � 30°a re given in Figures 17-20 . It is noticed that the impact of surface water infiltration on the required supporting pressure is more obvious with the increase of the angle of internal friction, with the maximum increment of 57 kPa when the friction angle is taken as 30°. e induced variation of critical tunnel spacing is, however, becoming smaller. e plastic dissipations of the tunnels in clay with a friction angle of 20°and the center-to-center spacing of 3D subject to different magnitudes of surface water infiltration are compared in Figure 21. It is observed that the interaction between the two tunnels appears under hydrostatic conditions, (d) (e) (f )      which is not shown after the surface water infiltrate into the surrounding soils. erefore, the effect of surface water infiltration and ground water table fluctuation have to be considered in design of dual tunnels to assure their longtime stability during construction and operation.

Conclusions
e required supporting pressures of dual circular tunnels in unsaturated soils are obtained in the present paper using a unified computational limit analysis method. For simplicity, a uniform pressure is exerted on the periphery of the tunnel to stabilize the surrounding soils. e hydraulic states in the cases of ground water table fluctuation and surface water infiltration are considered, and the effect of the pore water saturation of the surrounding soils on the loosening pressure, the collapse mechanism, and the critical center-tocenter spacing of dual tunnels is discussed. e following conclusions can be obtained as follows: (1) e established numerical algorithm is proved to be efficient and can be reliable for the stability analysis of dual tunnels in unsaturated soils. (2) For all considered ground water tables, the required supporting pressure decreases gradually with the increase of the tunnel spacing, before reaching a constant value. (3) e ground water table has a large impact on the required supporting pressure and the critical tunnel spacing of dual tunnels. Smaller pressures are needed to maintain the stability of the tunnels with the decrease of the ground water table, due to the additional shearing strength contributed by the matric suction. e effect of the ground water table is most significant when it is located near the tunnels. (4) e required supporting pressure increases if there is surface water infiltrating into the surrounding soils and the effect of water infiltration depends on the hydraulic properties and the friction angle of surrounding soils. e variation of the critical tunnel spacing induced by infiltration decreases with the increase of the friction angle.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.