Seismic Behavior of Triple Tunnel Complex in Soft Soil Subjected to Transverse Shaking

Combining multiple tunnels into a single tunnel complex while keeping the surrounding area compact is a complicated procedure. The condition becomes more complex when soft soil is present and the area is prone to seismic activity. Seismic vibrations produce sudden ground shaking, which causes a sharp decrease in the shear strength and bearing capacity of the soil. This results in larger ground displacements and deformation of structures located at the surface and within the soil mass. The deformations are more pronounced at shallower depths and near the ground surface. Tunnels located in that area are also affected and can undergo excessive distortions and uplift. The condition becomes worse if the tunnel area is larger, and, thus, the respective tunnel complex needs to be properly evaluated. In this research, a novel triple tunnel complex formed by combining three closely spaced tunnels is numerically analyzed using Plaxis 2D software under variable dynamic loadings. The effect of variations in lining thickness, the inner supporting structure, embedment depth on the produced ground displacements, tunnel deformations, resisting bending moments, and the developed thrusts are studied in detail. The triple tunnel complex is also compared with the rectangular and equivalent horizontal twin tunnel complexes in terms of generated thrusts and resisted seismic-induced bending moments. From the results, it is concluded that increased thickness of the lining, inner structure, and greater embedment depth results in decreased ground displacements, tunnel deformations, and increased resistance to seismic-induced bending moments. The comparison of shapes revealed that the triple tunnel complex has better resistance against moments with the least amount of thrust and surface heave produced.


Introduction
Tunnels are one of the important means of underground transportation. They are useful in the build-up area where further surface construction in order to accommodate large traffic infrastructure is not possible. Second, they also provide fast, uninterrupted flow from one point to the other. Many densely populated cities in the world are now planning to provide metro lines through an underground tunnel system. Although it has many useful purposes, at the same time, it is vulnerable to excessive damage if present in an earthquake-prone area. The condition becomes more critical when multiple tunnels lie in closer proximity. Previously, it was believed that the underground structures are less prone to damage in comparison with the on-surface constructions, but major earthquakes in the recent past like Kobe (1995), Coyote (1979), Kocaeli (1999), Chi-Chi (1999), etc. have proved these to be equally vulnerable. The damage suffered by the tunnels depends upon the type of surrounding soil, embedment depth, the amplitude of the earthquake, and the groundwater table (GWT) conditions. The type of soil plays a very important role in the overall behavior of tunnels under seismic activity. Tunnels constructed in soft soil are far vulnerable to severe damage than in dense rock. Similarly, the more shallow tunnels suffer more distortions in comparison to the deeper ones [1].
Researchers now focus on the dynamic behavior of tunnels in different types of rocks and soils. Wang [2], Penzien [3], and Bobet [4] developed closed-form solutions to estimate the seismicallyinduced moments and thrust in rectangular and circular tunnels causing racking and ovaling, respectively. Hashash and Hook et al. [5] and Hashash and Park et al. [6] conducted an extensive review of the available methods used to determine seismic induced forces that aid in the seismic design of different tunnels. Liu and Song [7] and Azadi and Hosseini [8] studied the seismic behavior of shallow tunnels subjected to horizontal and vertical shaking in the liquefiable sands in order to assess the distortions occurring in the tunnel lining and its uplift due to excess pore pressure generation. Unutmaz [9] investigated the liquefaction possibilities of soil surrounding the tunnels using 3D finite element modeling (FEM) and studied the acceleration variations and the power spectra developed under the seismic vibrations for various tunnel embedment depths. Lanzano and Bilotta et al. [10] studied the behavior of circular tunnels, both using a centrifuge and analytical modeling in order to determine the hoop forces and the distortions produced in the tunnel lining and concluded that both techniques produce results in good agreement. Cilinger and Madabhushi performed centrifuge tests to study the effect of embedment depth and behavior of circular and rectangular tunnels in dry sand subjected to different intensity seismic vibrations [11][12][13]. Qiu and Xie et al. [14] conducted centrifuge testing to study the effect of ground movements and the interaction between twin tunnels in loess and concluded on the optimized spacing and the interval between the tunnels. Shahrour and Khoshnoudian et al. [15] numerically studied the seismicallyinduced bending moments and soil dilatancy around tunnels in soft soils. Patil and Choudhury et al. [16] performed FEM on circular tunnels and studied the lining distortions, bending moments under variable ground vibrations, and different embedment depth and tunnel lining thickness to develop better understanding of tunnel behavior in soft soils. Tsinidis [17] conducted FEM to study response characteristics of rectangular tunnels under seismic shaking varying the tunnel-soil interface properties and input motion characteristics embedment depths and concluded in the development of racking-flexibility (RF) relations for rectangular tunnels in soft soils. Apart from this, the response of underground tunnels subjected to seismic ground shaking have been studied experimentally (Chian and Madabhushi [18], Graziani and Boldini [19], Chian and Madabhushi [20], Abuhajar and Naggar et al. [21]), analytically (Power et al. [1], Chian and Tokimatsu [22], Bobet and Fernandez et al. [23]) and numerically (Chou and Yang et al. [24], Amorosi and Boldini [25], Baziar and Moghadam et al. [26], Huo and Bobet et al. [27], Nguyen and Lee et al. [28]). However, all literature considers the conventional circular, rectangular or horseshoe shapes. None of them has studied the behavior of multiple tunnels that have been combined to form a single tunnel complex of an unconventional shape.
This research includes the dynamic response of a new, unconventional tunnel shape, which is proposed to be constructed in Brussels, which is the capital of Belgium, and has never been used earlier. The literature also lacks analytical solutions for it. This shape has resulted from the combination of three closely spaced tunnels into a triplet complex that would carry multiple train tracks for the underground traffic movement and has been discussed with respect to construction arrangements by Naseem and Schotte et al. [29]. This unconventional triplet complex in soft soil has been numerically studied in detail using the dynamic module of FEM software Plaxis 2D under variable seismic vibrations and embedment depths in order to determine the tunnel distortions, which results in ground settlements, seismically-induced bending moments and axial forces, and make the comparison with the conventional rectangular complex and the equivalent horizontal twin tunnel complex. This study would help determine the optimized tunnel parameters and to understand the dynamic behavior of this unique tunnel shape in a better way.

Stages of Analysis
The research is comprised of static and dynamic analyses. First, the tunnel system is analyzed statically. Since it is a long-term response, drained conditions are considered. After the equilibrium conditions are established, the dynamic analyses are performed, which include rapid ground shaking. Therefore, undrained conditions are considered at that stage.

Reference Model
This research focuses on the parametric study of an unconventional shape and a unique triple tunnel complex under seismic shaking in order to evaluate deformations in tunnel lining, moments, and thrust developed. As studied for the first time, it lacks the analytical solutions. Hence, it is compared with a reference model. The reference model considered in this case is taken from the research of Milind and Choudhury et al. [16]. It comprises of a circular tunnel with a diameter of 6 m analyzed using FEM software Plaxis 2D under different seismic vibrations. The Mohr-Coulomb constitutive model was used in the research and the results of the developed thrusts and moments were compared with the analytical solutions of Wang [2], Penizen [3], and Bobet [4] as shown in Figure 1. The results obtained were in close agreement except from one anomalous result of thrust in case of Penizen [3], which indicates their accuracy.
In this research, the seismic shaking of a triple tunnel complex is analyzed using the same FEM suite with the same layered soil system and the constitutive model so that the obtained results have the confidence of accuracy.

Software and Soil Parameters
In this study, the FEM program Plaxis 2D is used, which offers a wide range of constitutive models and different modules to numerically analyze both the linear and non-linear behavior of the soil and tunnel structure. A layered soil system is considered to have a soft silty clay layer underlain by silty clay, very soft silty clay, clay, and silty clay with silty sand, respectively, as used by Hu and Yue et al. [30] along with Milind and Choudhury et al. [16]. The soil media is taken as fully saturated with the GWT considered to be at the ground level. The triple tunnel complex is provided in the soil system keeping the embedment depth (C) and tunnel complex width (H) as a variable. The triple tunnel complex along with the layered soil system is shown by Figure 2a-c while the detailed soil properties used in the study are tabulated in Table 1. This layered soil system is categorized as soft soil with soil type D, according to Euro code (EC 8) [31] based on its shear wave profile given in Figure 2d.

Constitutive Model and Boundary Conditions
The soil tunnel system is modeled using Plaxis 2D, considering the Mohr-Coulomb constitutive model for the soils used by many other researchers [16,30,32]. In order to use for the dynamic analysis, Plaxis 2D makes use of the modified Mohr-Coulomb constitutive model. The shear, elastic moduli and other parameters in this model are dependent on the primary input of the shear wave velocity of the soil medium and, thus, are calculated accordingly. Apart from this, shear and elastic moduli variation w.r.t depth and confining stresses are also taken into account. The equations can be given by the equation below.

=
(1) where G is the shear modulus, Vs is the shear wave velocity, E is the elastic modulus, υs is the soil's Poisson's ratio, y is the depth in the vertical direction, Eref is the value of elastic modulus at reference point, yref is the depth of the reference point, and Einc is the increment in value of the elastic modulus [33,34]. The seismic ground shaking produces cyclic stresses in the soil, which generates a hysteresis loop with the dissipation of energy and damping. The Mohr-Coulomb soil model is an elasticperfectly plastic soil model and cannot capture this phenomenon. To cater this limitation, the where α and β are the Rayleigh viscous damping coefficients, ω is the angular frequency, h is the thickness of the soil layer, f1 and f2 are the first and second target frequencies, and ξ1 and ξ2 are the respective damping ratios. The damping ratio is taken as 10% for all of the soil layers that attributes to the soft soil [35]. The damping coefficients are calculated individually for each of the soil layers incorporated into the constitutive model definition. The tunnel lining and inner supporting structure are modeled using linear elastic plate elements as used by many other researchers in the past [8,9,16,17]. The Elastic modulus is taken as 37 MPa while the unit weight of concrete and the Poisson's ratio are taken as 25 kN/m 3 and 0.2, respectively. The lateral boundaries are kept far apart at a distance of 400 m based on the sensitivity analysis so that they do not interfere with the passage of earthquake waves. The tied degree of freedom condition is assigned to the lateral boundaries so that the left and right nodes move together while absorbing the lateral waves. This condition is used to simulate the freefield phenomenon of the soil. The base of the model is taken as rigid, which ensures that the waves do not escape and are reflected back into the model. The base is kept at a depth of 75 m and the earthquake signals are applied at the base of the model. The schematic diagram of the defined model geometry can be seen in Figure 2a. The model is discretized using 15-noded triangular elements with the mesh size selected based on Equation (7) given by Kuhlemeyer and Lysmer [36].
And λ = where Δl is the length of the finite elements and λ is the wavelength, f is the frequency, and Vs is the least shear wave velocity among all the soil layers. The time-step also plays an important role in the overall accuracy of the results and, hence, is selected by the equation below.
= * where δt is the time step, t is the total time duration of the seismic vibration. In addition, m and n are the maximum number of steps and the number of sub-steps, respectively. This ensures that the waves do not pass more than one element at each time step.

Input Ground Motion Characteristics
In order to evaluate the seismic performance of any structure, the major devastating earthquakes of higher magnitudes in history like Kobe (Japan, 1995), Loma Prieta (USA, 1989), Chi-Chi (Taiwan, 1999), and Coyote (USA, 1979) are used, which can be seen from the literature. In this study, two major earthquakes acceleration-time histories are used to analyze the dynamic behavior of the triple tunnel complex. The ground motions used are of the Coyote earthquake (USA, 1979) and the Kocaeli earthquake (Turkey, 1999). The acceleration-time histories, their Fourier response spectra, and spectral accelerations are shown in Figure 3a,b while their details are tabulated in Table 2. In order to use in the research, all of them are scaled to 0.4 g unless otherwise stated and named as input motion (IM) 1 and 2, respectively. The spectral accelerations of the input motions are also compared with the EC8 site class A and their average normalized curve compares well with the design spectrum, which can be seen in Figure 3c.

Parametric Study
This research is focused on the dynamic analysis of a novel triple tunnel complex, which requires an exhaustive study. In order to better understand its behavior, a detailed parametric study is conducted by varying the thickness of lining, inner supporting structures, and changing the embedment depth to a tunnel complex width (C/H) ratio under the effect of variable amplitudes of seismic vibrations to determine the effect of each of the parameters on the overall tunnel deformations and the resulting surface displacements. Apart from the parametric study, two other tunnel shapes named a rectangular complex and the equivalent horizontal twin complex are also analyzed for the seismic-induced bending moments and thrusts. The obtained results are compared with the triple tunnel complex to evaluate its performance. On the basis of the results obtained, the conclusions are then drawn.

Effect of the Variation in Lining Thickness
In order to study the effect of tunnel lining thickness, the tunnel lining thickness is kept variable from 0.4 m to 0.8 m, which constitutes 1.63% to 3.26% of the total tunnel width. All the other factors, i.e., inner structure's thickness, the amplitude of the IM and the C/H ratio are kept constant and the overall surface displacements are calculated. Apart from this, the normalized tunnel deformations and the bending moments are also computed. Normalized tunnel deformation is defined as the ratio of distortions between the top and bottom parts of the tunnel complex to the distortions at the same points in the free-field. This parameter is important in determining the flexibility of the tunnel lining. From the results shown in Figure 4b, it can be depicted that the normalized tunnel deformations are almost equal to 1 when the lining thickness is 1 m. It means that the lining thickness greater than 1 m would behave as rigid while a lesser thickness would act as flexible. The increased lining thickness causes a reduction in the flexibility of the tunnel and resistance to deformations increases. It also results in the reduced surface heave with increased resistance to seismic-induced bending moments,  [16] and by Azadi and Hosseini [8]. From the figures, it is evident that the tunnel complex with thicker lining has greater resistance to the seismic-induced uplift, distortions, and bending moments. From the study, it was also noted that reduction in lining thickness from 0.4 m to 0.3 m made the tunnel complex very fragile and resulted in the collapse of the structure.

Effect of Variation in the Inner Structure
The triple tunnel complex is formed by truncating some of the portions of the three closely spaced tunnels and joining them with the slab and the concrete walls in between, as shown in Figure  2c [29]. This inner structure plays an important role in the overall dynamic behavior of the triple tunnel complex and, hence, is studied in this research. The thickness is varied from 0.3 m to 0.7 m while keeping all the other parameters, i.e., lining thickness, C/H, and amplitude of the IM as constant. From Figure 5, it can be seen that, as the thickness of the inner supporting structure increases, the surface heave decreases. Even though the reduction is minimal in comparison to when the lining thickness is varied by the same amount, it can be inferred that the thicker inner supporting structure would add to the enhanced resistance of the triple tunnel complex. From the study, it was also noted that, when the thickness was further reduced to 0.2 m, the tunnel complex was unable to bear the seismic vibrations and resulted in the collapse.

Effect of the Variation in the Amplitude of the Input Motion
To study this effect, the IM 1 is scaled up and down in order to have a range varying from 0.1 g to 0.5 g. Keeping all other parameters, i.e., tunnel thickness, inner structure thickness, C/H as a constant, the IM is applied in order to compute the surface displacements and the corresponding tunnel uplift. From Figure 6a,b, it can be seen that, as the amplitude of the IM increases from 0.1 g to 0.5 g, the surface heave increases from 0.011 m to 0.105 m. The increase in the amplitude of IM results in a greater tunnel uplift as compared to the surface heave, which is minimal at the lower amplitudes but more pronounced at the higher amplitudes. This trend also follows a similar pattern as obtained by Azadi and Hosseini [8].

Effect of Variation in the Embedment Depth
To study the effect of embedment depth on the dynamic behavior of the triple tunnel complex, a tunnel with constant lining thickness and inner structure embedded at different C/H ratios was analyzed using a constant amplitude IM 2. Figure 7a, Ground surface Tunnel settlements and normalized tunnel deformations. From Figure 7a, it is evident that, as the embedment depth increases, the amount of surface heave decreases while the plot in Figure 7b shows that the increased embedment depth causes the normalized tunnel deformations to decrease. It means that the triple tunnel complex embedded deeper would have more resistance to the lining distortions than the shallower tunnel complex. In other words, a deeper embedded tunnel complex would behave as rigid in comparison to the same tunnel complex embedded shallower.

Effect of Variation in the Shape of the Tunnel Complex
To study the effect of shape on the dynamic behavior of the tunnel complex, two other shapes named rectangular and an equivalent horizontal twin tunnel complex [29], which are shown in Figure  8a,b, respectively, are also analyzed. Keeping all other parameters, i.e., tunnel lining, inner structure, embedment depth, and amplitude of the IM as constant, the thrusts and moments are computed for the previously mentioned shapes and a comparison is made. Figure 9a-c are showing the seismicinduced bending moments and thrusts in the triple tunnel complex, rectangular tunnel complex, and equivalent horizontal twin tunnel complex, respectively. In order to compare the dynamic performance, the surface displacements produced in the presence of each of the individual tunnel complexes and the variation of moments and thrusts are plotted along the normalized tunnel perimeter, which can be seen in Figure 10a-d. From Figure 10a, it is clear that the surface heave developed in the presence of a triple tunnel complex is the minimum among all, equaling 0.082 m, followed by 0.113 m by a rectangular tunnel complex while the maximum is produced by the equivalent horizontal twin tunnel complex, which is 0.135 m, i.e., almost 36.44% and 63.85% more, respectively. From Figure 10b,c, it can be noted that the resisted seismic-induced bending moments along the normalized perimeter of about 0.375, 0.5, and 0.625, the resisted moments in case of the triple tunnel complex are about 79.86% more than the rectangular tunnel complex. The maximum resisted moments compared to the equivalent horizontal twin complex are almost 6.55% more.

Conclusions
This research has been carried out to study the dynamic behavior of a triple tunnel complex in soft soil and the effect of different parameters on the overall behavior in case of transversal shaking by an earthquake. From the detailed study conducted, the following conclusions are drawn. the deeper embedded tunnel would act more rigid in comparison to the shallower embedded tunnel with the same lining thickness. (5). The shape of the tunnel complex also plays an important role in its seismic behavior. The triple tunnel complex results in a surface heave of 0.083 m as compared to 0.113 m and 0.136 m, which are 36.44% and 63.85% more in the presence of a rectangular tunnel complex and the equivalent horizontal twin tunnel complex, respectively. (6). The maximum resisted seismic-induced bending moments by the triple tunnel complex are 33.37% less than that of the rectangular tunnel complex while 6.55% more than the equivalent horizontal twin tunnel complex at the tunnel invert. Looking to the resisted seismic-induced bending moments along the tunnel perimeter, it can be noted that, at the normalized perimeter of about 0.375, 0.5, and 0.625, the resisted moments in case of the triple tunnel complex are about 30.20% more while, at the crown level, about 79.86% more than at the rectangular tunnel complex.
From the conclusions drawn, it is evident that the triple tunnel complex resists more seismicinduced bending moments than that of the equivalent horizontal twin tunnel complex and the rectangular tunnel complex. The produced surface heaves are the least among the three shapes as well. This highlights the performance of a triple tunnel complex against the transversal seismic vibrations.

Limitations
Since this study comprises of a new, unconventional shape that lacks analytical solutions, the reliability of results is based on the comparison of seismic-induced bending moments and thrusts results with the analytical solutions of the reference model. For more authenticity, analytical solutions addressing this particular shape are needed to develop. Apart from this, experimentation using centrifuge must be conducted in order to validate and compare the results. Currently, the literature lacks research regarding this unique shape. Hence, on the basis of the reference model, the Mohr-Coulomb model modified for incorporating the effect of cyclic stresses using viscous damping and moduli increment is used in order to compare the results. Since the earthquake also results in the shear modulus degradation with the development of shear strains, advanced models like the Hardening strain with small strain stiffness (HSsmall) or the soft soil model should also be used. The results should be compared and more refined.

Future Work
This research is focused on the detailed parametric study of the triplet tunnel complex and its comparison with the other tunnel shapes numerically. In this research, continuous tunnel lining is used with the minimal internal structure thickness for the analyses. Therefore, the conclusions are valid for that condition. It is recommended that the effect of segmental lining and fortification of the tunnel complex using a thicker internal structure should also be evaluated. Apart from this, the response spectrum analysis should be performed to further evaluate its performance.