#### Abstract

The 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. The linings are substituted by an equivalent uniform pressure exerted on the periphery of the tunnel. The 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. The critical center-to-center spacing above which the interaction of the tunnels disappears is also discussed.

#### 1. Introduction

The 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. The 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. The 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.

The ground movements around parallel tunnels in clays were investigated by Wu and Lee [1] by a series of centrifuge model tests. The 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 earth-pressure-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–9]. Recently, CLA was frequently applied in the stability analysis of twin tunnels under the condition of plain strain [10–16], and the emphasis was mainly given to the impact of the factors such as critical tunnel spacing and novel section shapes. The 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. The 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. The matric suction is affected by the external hydraulic conditions such as rainfall, irrigation, and ground water level fluctuation. Therefore, 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–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. The 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.

The 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-to-center spacings of tunnels in three different soils. Concluding remarks are provided in Section 6.

#### 2. 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.

##### 2.2. Upper Bound Theory for Unsaturated Soils

CLA is a direct method for determination of limit loads based on the classical bounding theory. The 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, is the pore air pressure, is an identity vector, and is the suction stress. The closed form solution of the suction stress is [24]where is the pore water pressure and and 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. The soils are then dealt with as equivalent single-phased materials and the strength parameters are invariant.

The following mathematical programming problem is established to obtain the limit loads of unsaturated earth structures:where is the load multiplier, is the strain rate, is the velocity, and are the unified body force and the unified surface force specified beforehand, 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

It is observed that equation (4) turns into the formulation for limit analysis of saturated soils when pore water pressure is considered. Therefore, unsaturated soils could be consistently analyzed by the present model.

##### 2.3. Numerical Discretization

The 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. The velocities in an element are expressed in terms of the element velocity vector aswhere and are the components of the velocity and is a coefficient matrix consisting of shape functions of the velocities. Correspondingly, the strain rate in some element is given aswhere is the area coordinate with respect to the *i*th vertex, is a coefficient matrix consisting of shape functions of the strains, and is the element strain rate vector. The effective stress within an element iswhere is the coefficient matrix for the effective stresses, is the element effective stress vector, and is the effective stress at the vertexes. The strain rate within an element can be obtained by the element velocity vector aswhere is the strain matrix and is a differential operator relating the strain rate and the velocities. Introduction of equations (5)–(8) into equation (3) yieldswhere is the number of elements and is the number of element boundaries where surface tractions are specified. Equation (9) can be simplified aswhere is the global effective stress vector, is the global external force vector, and 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. The 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.

##### 2.4. Second-Order Cone Programming

The soils are supposed to satisfy the Mohr–Coulomb yield criterion:where and are the tangential and normal effective stresses on the failure surface and and are the effective cohesion and angle of internal friction, respectively. The minus sign in front of the normal stress appears because the stress is supposed to be negative in compression. The 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 is firstly introducedwhere is the mean effective stress, and are the components of the effective stress, is the deviatoric stress, and is Kronecker’s delta. For the condition of plain strain, the failure criterion is rewritten as

An auxiliary variable *z* is introduced to reformulate equation (13) intowhere and . The discretized unified upper bound formulation will then be summarized aswhere 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.

#### 3. 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. The 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 and 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 and the matric suction:where and have the same meanings as those in equation (2).

#### 4. Verification

The 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. The 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. The obtained normalized bearing capacities are compared with those in [21] in Figure 2. The 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.

#### 5. Stability of Dual Tunnels

##### 5.1. Problem Definition

The configuration of the dual tunnels studied in the present paper is given in Figure 3. The 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. The buried depth and the diameter of the tunnels are *H* and *D*, and the center-to-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. The 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. This 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. The distance between the ground water table and the top ground is . The water table will be placed at some typical locations, indicated as *P*_{A} − *P*_{J}, as shown in Figure 3(a). The points of *P*_{H}, *P*_{I}, and *P*_{J} are not given in the figure, and the vertical intervals of those points are *D*. The selected locations include the top ground, the top of the tunnel, the middle of the tunnel, the bottom of the tunnel, and other positions with the vertical distance of *D* or *D*/2. For the symmetry of the problem, half of the solution domain is considered, as shown in Figure 3(b). Horizontal displacement of the left boundary and all the displacements of the bottom and right side of the solution domain are fixed. A standard mesh is given in Figure 3(c), where it is observed that a large enough solution domain is chosen to eliminate the influence of boundary conditions and a finer mesh is utilized near the tunnel to increase computational efficiency. Before parameter analysis, the authors have carried out some preliminary analyses to check the size of the solution domain and the chosen of the element mesh. It is shown that the results are not affected by the solution domain and the element mesh.

**(a)**

**(b)**

**(c)**

##### 5.2. 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. The dimensionless stability number is defined as *N* = _{lim}*D*/*c*, where _{lim} is the maximum unit weight of the surrounding soils the tunnels can resist without support. The 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.

##### 5.3. Effect of Ground Water 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. Therefore, 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. The largest suction stress is seen in clay, whereas the suction stress in sand is the smallest. The 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. The 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.

In the following discussions, hydraulic parameters of the surrounding soils are taken as those of silt. The effective cohesion is 10 kPa and the effective angle of internal friction is chosen to be 10°, 20°, and 30°, respectively. The 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. There 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 5*D* to 3.5*D* 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. The 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. The 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 2*D* and the ground water table is placed at *P*_{G}. The negative value of the limit supporting pressure means the tunnel can maintain its stability without any support measures. Therefore, 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.

**(a)**

**(b)**

**(c)**

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. The 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 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. The situation is different in the case of *S*/*D* = 3, where it can be noticed that, as /*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 /*D* is larger than 2. Therefore, 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.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

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. The same form of variation is observed as that in the case of *H*/*D* = 1. The critical center-to-center tunnel spacing increases with the increase of the buried depth of the tunnels. The 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. Therefore, the buried depth has to be considered in determination of the critical tunnel spacing. The 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. The 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. Therefore, it can be appropriate to say that the drawdown of the ground water table aggravates the influence of tunnel spacing.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

##### 5.4. 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. The 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. The 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]where is the infiltration rate and is the unit weight of the pore water. The 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. The extent of reduction is related to the type of the soil. The 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. The 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. The 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.

**(a)**

**(b)**

**(c)**

Variations of the limit supporting pressures for = 10° are 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. The critical center-to-center spacing is reduced by the water infiltrating into the surrounding soils. The 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. The limit supporting pressure grows fast at the instant of infiltration, and then, the changing rate almost reaches a constant value. The dual tunnels in different surrounding soils are affected to different degrees by water infiltration. The maximum difference of critical spacing caused by infiltration 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.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

The limit supporting pressures for = 20° and = 30° are 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°. The induced variation of critical tunnel spacing is, however, becoming smaller. The plastic dissipations of the tunnels in clay with a friction angle of 20° and the center-to-center spacing of 3*D* 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, which is not shown after the surface water infiltrate into the surrounding soils. Therefore, the effect of surface water infiltration and ground water table fluctuation have to be considered in design of dual tunnels to assure their long-time stability during construction and operation.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

#### 6. Conclusions

The 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. The 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-to-center spacing of dual tunnels is discussed. The following conclusions can be obtained as follows:(1)The 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)The 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. The effect of the ground water table is most significant when it is located near the tunnels.(4)The 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. The variation of the critical tunnel spacing induced by infiltration decreases with the increase of the friction angle.

#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was performed with the support of Fundamental Research Funds for the Central Universities (no. 300102210209) and Natural Science Basic Research Plan in Shaanxi Province of China (no. 2018JQ5098).