Analysis of Undrained Seismic Behavior of Shallow Tunnels in Soft Clay Using Nonlinear Kinematic Hardening Model

: In this study, a soil–tunnel model for clay under earthquake loading is analyzed, using ﬁnite element methods and a kinematic hardening model with the Von Mises failure criterion. The results are compared with those from the linear elastic–perfectly plastic Mohr–Coulomb model. The latter model does not consider the sti ﬀ ness degradation caused by imposing cyclic loading and unloading to the soil, whereas the kinematic hardening model can simulate this sti ﬀ ness degradation. The parameters of the kinematic hardening model are calibrated based on the results of experimental cyclic tests and ﬁnite element simulation. Here, two methods—one using data from cyclic shear tests, and the other a new method using undrained cyclic triaxial tests—are used to calibrate the parameters. The parameters investigated are the peak ground acceleration (PGA), tunnel lining thickness, tunnel shape, and tunnel embedment depth, all of which have an e ﬀ ect on the resistance of the shallow tunnel to the stresses and deformations caused by the surrounding clay soils. The results show that unlike traditional models, the nonlinear kinematic hardening model can predict the response reasonably well, and it is able to create the hysteresis loops and consider the soil sti ﬀ ness degradation under the seismic loads. Zhang’s equations [32] are used. The values of the parameters such as the undrained shear strength ( 𝑆 𝑢 ), maximum shear modulus ( 𝐺 𝑚𝑎𝑥 ) or shear wave velocity ( 𝑉 𝑠 ), and curve of the undrained Young’s modulus versus axial strain ( 𝐸 𝑢 Vs. 𝜀 1 ) are required to calibrate the parameters. A series of undrained cyclic triaxial and monotonic tests of kaolin clay were conducted by Wichtmann [33] in both stress and strain control modes and under different level stress amplitudes. this study, results of monotonic and cyclic tests under the stress control condition mean of kPa


Introduction
With the growing populations in metropolitan cities around the world, underground urban infrastructure development has become a major focus. If these underground structures are constructed in earthquake areas, their dynamic response to seismic loads must be determined in order to reduce economic damage and threats to human safety. In contrast to overground structures where design is based on the fact that earthquakes manifest as inertial forces, in underground structures such as tunnels, which have a significant degree of containment, design and analysis are based on the assumption that earthquakes primarily deform the soil structure. Previously, engineers believed that earthquakes caused less destruction of underground structures than surface structures, but the earthquakes of the 1990s caused serious damage to tunnels, forcing a re-appraisal of this view. In recent years, there have been numerous examples in different countries of earthquakes causing extensive damage to underground structures (e.g., the Kobe earthquake, 1995; Taiwan, 1999) [1][2][3][4][5][6][7], which indicates that it is crucial for studies to be made of the seismic behavior of tunnels.
As earthquake waves are usually amplified by passing through different layers of soil (e.g., clay or sand) [8], the surrounding soil can affect the performance of structures such as tunnels [9,10]. Previous studies have shown that tunnels located in clay soils are more likely to be damaged by earthquakes than those constructed in sandy soils [11][12][13][14]. However, there are fewer studies on the seismic behavior of tunnels located in soft soils than in sandy soils. Different approaches such as numerical [15][16][17] and experimental methods [18,19] are used to investigate the seismic behavior of tunnels in terms of the interactions between soil and tunnel.
Using numerical methods to investigate the dynamic behavior of tunnels requires a suitable constitutive model for the materials involved. In recent years, researchers have used numerical methods with different constitutive models to study the seismic behavior of tunnels. Kontoe et al. [5] used a modified Cam-clay constitutive model to investigate the dynamic response of tunnels during the 1999 Duzce earthquake in Turkey. Their results showed that this constitutive model was not suitable for simulating the seismic behavior of a tunnel during strong earthquakes. Debiasi et al. [20] investigated the dynamic response of shallow rectangular tunnels to earthquake loads. Their analysis was conducted mainly with pseudo-static finite element analyses, and it showed that these structures are strongly affected by nonlinear frictional effects at the soil-structure interface. Patil et al. [21] conducted a numerical study on the behavior of tunnels under seismic loads using PLAXIS2D software and a modified Cam-Clay model to model the soil medium. The results of the nonlinear plane strain analysis of the soil-tunnel system showed that the earthquake frequency affected the maximum and residual dynamic Earth pressure around the tunnel. Lei Zhang and Yong Liu [13] studied the seismic behavior of tunnel in clays under earthquake loads, with consideration of the hysteresis-hyperbolic constitutive model for soil. Their results suggested that randomness in the small-strain shear modulus of clay can significantly influence the shear force, the seismic responses of bending moments, and the lateral deformation of tunnels. Xiaohua Bao et al. [22] studied the dynamic behavior of a metro subway tunnel under earthquake loads, using a cyclic elastoplastic constitutive model for foundation soils so as to obtain acceptable results, and they also studied various parameters that affect the seismic behavior of a tunnel. Their work showed that the acceleration response decreased at different depths in the model due to the phenomenon of liquefaction. Moreover, the increase in the amplitude of vertical acceleration increased the liquefaction area.
Soil response to loads is strain-dependent. The behavior of soil under very small strains (γ < 10 −5 ) is linear elastic, but under higher strains, its behavior is nonlinear. Therefore, it is essential to select an appropriate constitutive model to consider nonlinear strain-dependent stiffness under different loads in a numerical analysis. Although many scholars have used numerical simulation methods to investigate the seismic behavior of soil-structure interactions (e.g., soil-tunnel interactions), an appropriate and more advanced constructive model would better predict the actual behavior of soil under earthquake forces. The use of some constitutive models such as Prevost's multi-yield model and the intergranular strain anisotropy model (ISA) requires complex calibration to determine their various parameters, which decreases their utility in research [23]. However, isotropic/kinematic plastic hardening models, comprising isotropic and kinematic hardening components, accurately predict soil behavior under dynamic loads.
The main theories on the hardening of soil are isotropic hardening, the linear and nonlinear kinematic hardening model of Armstrong-Frederick [24], and the nonlinear kinematic hardening model of Chaboche [25]. Isotropic hardening occurs when a yield surface is uniformly large, whereas kinematic hardening occurs when the center of the yield surface moves in stress space without the initial yield surface size changing, which is known as the Bauschinger effect [26,27]. Finally, if the size and center of the yield surface change simultaneously, this is denoted as combined hardening. As soils exhibit anisotropic kinematic hardening behavior, isotropic models are not suitable for studying soil behavior under seismic loads. Thus, we use a nonlinear kinematic hardening model in this study, which is a simplified version of the combined isotropic/kinematic hardening model generated by assuming that the size of the yield surface is constant [28].
In this study, a kinematic hardening constitutive model based on the Von Mises failure criterion is been used to investigate the seismic behavior of tunnels in clay soil. In some constitutive models, such as the Cam-clay and modified Cam-clay models used in previous studies, the yield surface is presumed to expand with increasing loading level. This means that when the soil experiences a larger stress surface under dynamic loading, the behavior of the materials will be quite elastic at lower loading Appl. Sci. 2020, 10, 2834 3 of 22 levels until the soil sample again experiences stresses greater than the previous yield stress. Therefore, these constitutive models cannot model the Bauschinger effect nor determine the hysteresis loops that exist under cyclic loading, which means that these models cannot predict the actual behavior of soil under kinematic loads, such as seismic loads. However, a kinematic hardening model can be used as an advanced model for the simulation of soil behavior, thereby overcoming these problems. Here, the phenomenon of liquefaction has not been considered and there is no generation of pore pressures. The phenomenon of cyclical degradation and degradation soil stiffness are due to soil structure.
We investigate the factors affecting the behavior of tunnels in clay under earthquake load, including the earthquake intensity, tunnel shape, embedment depth, and lining thickness. Then, we compare our results with the results of Patila et al. [29]. The earthquake used in this study is a 1D earthquake with a dominant frequency of 1 Hz, and the 0.4 g and 0.83 g maximum accelerations that occurred in 1995 in Kobe, Japan. A viscous spring boundary is used around the soil model to prevent the return of earthquake waves into the model, thus ensuring that actual results are obtained.

Constitutive Model
The function of yield surface, considering the Von-Mises criterion, is defined by Equation (1) [30]: where σ o is the initial size of the yield surface under zero plastic strain, i.e., at the elastic stress limit, α ij is the back stress that refers to the center position of the yield surface, and J 2 σ ij is the second deviatoric stress invariant, given by Equation (2): where σ ij and α ij are the stress and back stress tensor, respectively, and S mn and t mn are the stress and back stress deviatoric tensor, respectively. At an infinitely large plastic strain, when σ ij approaches the effective maximum yield stress (σ y ) (Equation (3)), the back stress tensor (α ij ) is given by Equation (4), as follows: The important aspect of the kinematic hardening model is the back stress evolution law, which describes the translation of the yield surface in the stress space. The kinematic hardening model consists of two parts-a purely kinematic term and a relaxation term-and can be written as (Equation (5)): where C and γ o are material parameters, dε pl is the plastic strain rate tensor, and dε pl is the equivalent plastic strain rate. In relaxation terms, the back stress is related to the equivalent plastic strain and the linear kinematic model is transformed to the nonlinear kinematic model. The plastic flow rule associated with the yield function is as follows (Equation (6)) [30]: where the constant dλ is defined as dλ = dε pl . The detailed expression can be found in Appendix A. Parameter σ 0 represents the initial size of the yield surface at zero plastic strain and controls the onset of nonlinear soil behavior; it is defined as the fraction (λ) of the (σ y ) effective maximum yield stress at large plastic strain. Decreasing this ratio improves the model nonlinearity, but if it is too small, the model will be non-economic. In contrast, an excessive increase in this quantity has a negative effect on the model that leads to unrealistic results. In the 1D condition, it is given by Equation (7): where λ is the material constant. The evolution Equation of kinematic hardening can be integrated to give Equation (8): where α is the back stress, ε pl is the equivalent plastic strain, C is the initial kinematic hardening modulus (C = E = 2(1 + ν)G 0 ) that defines the maximum translation of the yield surface, and γ o defines the rate decrease of the kinematic hardening modulus. Parameter C can be derived by differentiating Equation (8) to give Equation (9): The occurrence of displacement of the yield surface and the creation of back stress (α ij ) in cyclic loading can be understood by considering the kinematic hardening model. At infinite large plastic strain (when σ ij approaches σ y ), the back stress (α ij ) is equal to α s = C γ o and (σ ij − α ij ) tends to σ 0 (see Figure 1), which means that Equation (5) tends to zero. 1 illustrates the parameters and notation used in the equations.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 23 stress at large plastic strain. Decreasing this ratio improves the model nonlinearity, but if it is too small, the model will be non-economic. In contrast, an excessive increase in this quantity has a negative effect on the model that leads to unrealistic results. In the 1D condition, it is given by Equation (7): where is the material constant. The evolution equation of kinematic hardening can be integrated to give Equation (8): where is the back stress, is the equivalent plastic strain, is the initial kinematic hardening modulus ( = = 2(1 + ) 0 ) that defines the maximum translation of the yield surface, and defines the rate decrease of the kinematic hardening modulus. Parameter can be derived by differentiating Equation (8) to give Equation (9): The occurrence of displacement of the yield surface and the creation of back stress ( ) in cyclic loading can be understood by considering the kinematic hardening model. At infinite large plastic strain (when approaches ), the back stress ( ) is equal to = ⁄ and ( − ) tends to 0 (see Figure 1), which means that Equation (5) tends to zero. 1 illustrates the parameters and notation used in the equations.

Calibration of Parameters and Model Validation
The parameters of constitutive models of soils are often obtained by comparing experimental results and numerical simulations using finite elements at the element test level [31]. In this study, experimental data from an undrained cyclic triaxial test are used to calibrate the parameter , and two methods for calibrating the parameters , , and 0 are presented. First, the curve of the undrained Young's modulus versus axial strain is drawn using the results obtained from the cyclic triaxial test and numerical simulations of the test to obtain the parameter. Second, data from a cyclic shear test and Ishibashi and Zhang's equations [32] are used.
The values of the parameters such as the undrained shear strength ( ), maximum shear modulus ( ) or shear wave velocity ( ), and curve of the undrained Young's modulus versus axial strain ( Vs. 1 ) are required to calibrate the parameters. A series of undrained cyclic triaxial and monotonic tests of kaolin clay were conducted by Wichtmann [33] in both stress and strain control modes and under different level stress amplitudes.
In this study, the results of monotonic and cyclic tests under the stress control condition with an effective mean pressure of 200 kPa are used to calibrate the constitutive model parameters. The process mentioned above, by which the parameters of the model are calibrated, is fully described below.

Calibration of Parameters and Model Validation
The parameters of constitutive models of soils are often obtained by comparing experimental results and numerical simulations using finite elements at the element test level [31]. In this study, experimental data from an undrained cyclic triaxial test are used to calibrate the parameter C, and two methods for calibrating the parameters λ, γ o , and σ 0 are presented. First, the curve of the undrained Young's modulus versus axial strain is drawn using the results obtained from the cyclic triaxial test and numerical simulations of the test to obtain the λ parameter. Second, data from a cyclic shear test and Ishibashi and Zhang's equations [32] are used.
The values of the parameters such as the undrained shear strength (S u ), maximum shear modulus (G max ) or shear wave velocity (V s ), and curve of the undrained Young's modulus versus axial strain (E u Vs. ε 1 ) are required to calibrate the parameters. A series of undrained cyclic triaxial and monotonic tests of kaolin clay were conducted by Wichtmann [33] in both stress and strain control modes and under different level stress amplitudes.
In this study, the results of monotonic and cyclic tests under the stress control condition with an effective mean pressure of 200 kPa are used to calibrate the constitutive model parameters. The process mentioned above, by which the parameters of the model are calibrated, is fully described below.

Determination of C
We use the results of the undrained cyclic triaxial tests at different stress amplitudes (C1-C8 and C32-C36) [33] to determine the nonlinear undrained Young's modulus. For instance, Figure 2 shows the deviatoric stress (q) and undrained Young's modulus (E u ) versus the axial strain (ε 1 ) curve with stress amplitudes of 60 kPa and 70 kPa for undrained cyclic triaxial tests of C7 and C8, respectively. Each loading and unloading loop (Figure 2a Figure 2 depicts the stiffness degradation of soil under cyclic loading and unloading. As shown in the diagrams, by increasing the number of cycles and axial strain amplitudes, the undrained Young modulus or the shear modulus decrease, and the amount of damping caused by the soil plasticity increases. Although by increasing number of cycles the amount of axial strain is also increased, the increase in axial strain on both sides of the horizontal coordinate axis is unequal, such that there is more on the positive side (compression). We use the results of the undrained cyclic triaxial tests at different stress amplitudes (C1-C8 and C32-C36) [33] to determine the nonlinear undrained Young's modulus. For instance, Figure 2 shows the deviatoric stress ( ) and undrained Young's modulus ( ) versus the axial strain ( 1 ) curve with stress amplitudes of 60 kPa and 70 kPa for undrained cyclic triaxial tests of C7 and C8, respectively. Each loading and unloading loop (Figure 2a,c) is used to calculate the undrained Young's modulus value at the desired strain amplitude, using the equation = 1 (Figure 2b,d). Therefore, each point in Figure 2b,d is the value of the undrained Young's modulus of the corresponding cycle in Figure 2a,c. Thus, Figure 2 depicts the stiffness degradation of soil under cyclic loading and unloading. As shown in the diagrams, by increasing the number of cycles and axial strain amplitudes, the undrained Young modulus or the shear modulus decrease, and the amount of damping caused by the soil plasticity increases. Although by increasing number of cycles the amount of axial strain is also increased, the increase in axial strain on both sides of the horizontal coordinate axis is unequal, such that there is more on the positive side (compression). ; (c,d) test C8 with = kPa , = kPa , = .
[33]. ( 1 ) based on the results of numerical analysis, Equation (11), and the experimental results of the several undrained cyclic triaxial tests, where these are performed under different stress amplitudes (C1-C8 and C32-C36) [33] to calibrate the parameters and of the constitutive model. The calibration of parameter using the data in Figure 3 is described in Section 2.2.2. Each experiment test [33] is based on a number of specified cycles, and numerical simulations were performed based on the same number of cycles. The parameters of the model have been calibrated in order to best fit the numerical model points with the experimental results. In this figure, the second and fifth cycles of curves of deviatoric stress ( ) versus axial strain ( 1 ) of the undrained cyclic triaxial tests are used  based on the results of numerical analysis, Equation (11), and the experimental results of the several undrained cyclic triaxial tests, where these are performed under different stress amplitudes (C1-C8 and C32-C36) [33] to calibrate the parameters C and λ of the constitutive model. The calibration of parameter λ using the data in Figure 3 is described in Section 2.2.2. Each experiment test [33] is based on a number of specified cycles, and numerical simulations were performed based on the same number of cycles. The parameters of the model have been calibrated in order to best fit the numerical model points with the experimental results. In this figure, the second and fifth cycles of curves of deviatoric stress (q) versus axial strain (ε 1 ) of the undrained cyclic triaxial tests are used to obtain the undrained Young's modulus values (E u ); then, these values are divided by the function of the void ratio (Equation (10)) and plotted versus the corresponding strain amplitude ε ampl 1 .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 23 to obtain the undrained Young's modulus values ( ); then, these values are divided by the function of the void ratio (Equation (10)) and plotted versus the corresponding strain amplitude 1 . The equation of the void ratio (used in Figure 3) obtained on the basis of experimental data is presented as follows [33] (Equation (10)): where is the void ratio of soil. As shown in Figure 3, by increasing the strain amplitude, the values of the undrained Young's modulus decrease, and this decreasing trend is proportional to Equation (11) [34]: where 1, = 3 * 10 −4 , corresponding to = , = 0.5 and = 0.9 (for kaolin clay) [33].
In addition, by subjecting the experimental data to extrapolation operations, e.g., ≈ 30 * , , Equations (9)-(11) can be used to express parameter as follows (Equation (12)): Using the results obtained from the undrained cyclic triaxial test at very small strains ( 1 < 10 −6 ), the value of the maximum undrained Young's modulus can be calculated.
Then, according to the equation , the value of the maximum shear modulus can be obtained, where and are the maximum undrained Young's modulus and Poisson ratio, respectively, and = 0.49 in the undrained case.

Determination of , σ 0 , and γ o
The undrained shear strength ( ) is obtained from the results of a monotonic consolidated undrained triaxial test. Figure 4 shows the stress-strain relationship in the monotonic undrained triaxial test [33]. According to this figure and using the following equation, the maximum value of deviatoric stress and undrained shear strength can be calculated as follows (Equation (13)). The Equation of the void ratio (used in Figure 3) obtained on the basis of experimental data is presented as follows [33] (Equation (10)): where e is the void ratio of soil. As shown in Figure 3, by increasing the strain amplitude, the values of the undrained Young's modulus decrease, and this decreasing trend is proportional to Equation (11) [34]: where ε ampl 1,r = 3 * 10 −4 , corresponding to E u E umax = f E u ,ampl = 0.5 and α = 0.9 (for kaolin clay) [33]. In addition, by subjecting the experimental data to extrapolation operations, e.g., E u max ≈ 30 * f E u ,e , Equations (9)-(11) can be used to express parameter C as follows (Equation (12)): Using the results obtained from the undrained cyclic triaxial test at very small strains (ε ampl 1 < 10 −6 ), the value of the maximum undrained Young's modulus can be calculated.
Then, according to the Equation G max = E umax 2(1+υ max ) , the value of the maximum shear modulus can be obtained, where E u max and υ max are the maximum undrained Young's modulus and Poisson ratio, respectively, and υ max = 0.49 in the undrained case. The undrained shear strength (S u ) is obtained from the results of a monotonic consolidated undrained triaxial test. Figure 4 shows the stress-strain relationship in the monotonic undrained triaxial test [33]. According to this figure and using the following equation, the maximum value of deviatoric stress and undrained shear strength can be calculated as follows (Equation (13)).
where S u is the undrained shear strength, σ 1 and σ 3 are the maximum and minimum mean stress, respectively, and the deviatoric stress is q = σ 1 − σ 3 .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 23 where is the undrained shear strength, 1 and 3 are the maximum and minimum mean stress, respectively, and the deviatoric stress is = 1 − 3 . The effective maximum yield stress considering the triaxial test condition is calculated as follows (the detailed expression can be found in Appendix B): The effective maximum yield stress considering the pure shear test is calculated as follows: The next parameter for calibration is . The first method of calibrating in the kinematic hardening model is based on the undrained cyclic triaxial test data. The deviatoric stress ( ) versus axial strain ( 1 ) under the stress control condition at different stress amplitudes are plotted based on the simulation results of the undrained cyclic triaxial test, and similar to the experimental method, the undrained elastic modulus ( ) of the second and fifth cycles of each test is obtained. Then, the curve of the undrained Young's modulus versus axial strain is plotted for both the experimental data and simulation results to calibrate . The simulation results of the undrained cyclic triaxial test and their comparison with the experimental results are presented in Figure 3. Here, = 0.14 is found to provide a reasonable fit to the experimental results, where the simulation is performed using Abaqus software under the stress control condition.
If undrained cyclic triaxial test data are not available, the results of the cyclic shear test under the strain control condition and the Ishibashi and Zhang experimental equations using the plasticity index of soil are used to calibrate the parameter λ (second method). In this method, the curve data of shear modulus versus shear strain ( − curve) is required to calibrate this parameter. Thus, the hysteresis loops of the cyclic shear test simulation in different strain amplitudes should first be plotted. For instance, the hysteresis loops of three different strain levels are shown in Figure 5. As can be seen, the value of the damping caused by the soil plasticity increases, and the shear modulus decreases if the strain levels are increased. Then, the shear modulus corresponding to the stabilized cycle is obtained for each strain amplitude to enable the plotting of the − points to calibrate (see Figure 6). For example, the shear stress-strain loops of three strain amplitudes in Figure 5 correspond to the designated points in Figure 6. These two figures illustrate the validation of the kinematic hardening model against the published − curve of Ishibashi and Zhang [32]. After The effective maximum yield stress considering the triaxial test condition is calculated as follows (the detailed expression can be found in Appendix B): The effective maximum yield stress considering the pure shear test is calculated as follows: The next parameter for calibration is λ. The first method of calibrating λ in the kinematic hardening model is based on the undrained cyclic triaxial test data. The deviatoric stress (q) versus axial strain (ε 1 ) under the stress control condition at different stress amplitudes are plotted based on the simulation results of the undrained cyclic triaxial test, and similar to the experimental method, the undrained elastic modulus (E u ) of the second and fifth cycles of each test is obtained. Then, the curve of the undrained Young's modulus versus axial strain is plotted for both the experimental data and simulation results to calibrate λ. The simulation results of the undrained cyclic triaxial test and their comparison with the experimental results are presented in Figure 3. Here, λ = 0.14 is found to provide a reasonable fit to the experimental results, where the simulation is performed using Abaqus software under the stress control condition.
If undrained cyclic triaxial test data are not available, the results of the cyclic shear test under the strain control condition and the Ishibashi and Zhang experimental equations using the plasticity index of soil are used to calibrate the parameter λ (second method). In this method, the curve data of shear modulus versus shear strain (G − γ curve) is required to calibrate this parameter. Thus, the hysteresis loops of the cyclic shear test simulation in different strain amplitudes should first be plotted. For instance, the hysteresis loops of three different strain levels are shown in Figure 5. As can be seen, the value of the damping caused by the soil plasticity increases, and the shear modulus decreases if the strain levels are increased. Then, the shear modulus corresponding to the stabilized cycle is obtained for each strain amplitude to enable the plotting of the G − γ points to calibrate λ (see Figure 6). For example, the shear stress-strain loops of three strain amplitudes in Figure 5 correspond to the designated points in Figure 6. These two figures illustrate the validation of the kinematic hardening model against the published G − γ curve of Ishibashi and Zhang [32]. After calibration using this method, λ = 0.12 is obtained.  After determining the parameters , , and , the initial size of yield surface ( 0 ) is obtained from Equation (7) ( 0 = ). The last parameter to be determined is the rate decrease of the kinematic hardening modulus ( ). According to Figure 1, the amount of effective maximum yield stress is determined by Equation (16): Finally, from Equations (14)- (16), parameter is calculated according to Equations (17) and  After determining the parameters , , and , the initial size of yield surface ( 0 ) is obtained from Equation (7) ( 0 = ). The last parameter to be determined is the rate decrease of the kinematic hardening modulus ( ). According to Figure 1, the amount of effective maximum yield stress is determined by Equation (16): After determining the parameters S u , σ y , and λ, the initial size of yield surface (σ 0 ) is obtained from Equation (7) (σ 0 = λσ y ). The last parameter to be determined is the rate decrease of the kinematic hardening modulus (γ o ). According to Figure 1, the amount of effective maximum yield stress is determined by Equation (16): Finally, from Equations (14)- (16), parameter γ o is calculated according to Equations (17) and (18): (Triaxial test condition) (Simple pure shear test condition). Generally, the triaxial test is an acceptable method for assessing the cyclic behavior of soils, and usually this test is used to calibrate the parameters of different constitutive models. Applying relatively simple boundary conditions and creating proper meshes in simulations of this test afford logical and reliable results. Thus, the results of the cyclic triaxial test using nonlinear kinematic hardening and Mohr-Coulomb models on a clay sample are compared with the experimental model results to validate and investigate the kinematic hardening constitutive model. This test is modeled in the 2D axisymmetric case to achieve acceptable results. Application of a 200 kPa confining pressure is simulated as a uniformly distributed load applied normally inward to all external elements around the specimen; then, the cyclic load (q = 60 kPa & 70 kPa) is applied to the top of the sample.
In the Mohr-Coulomb model, since the test is carried out in an undrained condition (i.e., an undrained cyclic triaxial test), the internal friction angle of soil (ϕ) is zero and the cohesion value (C) equals the undrained shear strength of soil (S u ). Figure 7 shows the comparison of deviatoric stress versus axial strain in numerical simulations (kinematic hardening model and Mohr-Coulomb) with the experimental results in stress control mode. As shown in the diagrams, the behavior of undrained soils using the kinematic hardening model is in relatively good agreement with the results obtained from the experimental test, while the Mohr-Coulomb model fails to create the hysteresis loops, and the results thus differ from the experimental results. In addition, the material behavior in the Mohr-Coulomb model is in the elastic range.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 23 (Triaxial test condition) (Simple pure shear test condition). Generally, the triaxial test is an acceptable method for assessing the cyclic behavior of soils, and usually this test is used to calibrate the parameters of different constitutive models. Applying relatively simple boundary conditions and creating proper meshes in simulations of this test afford logical and reliable results. Thus, the results of the cyclic triaxial test using nonlinear kinematic hardening and Mohr-Coulomb models on a clay sample are compared with the experimental model results to validate and investigate the kinematic hardening constitutive model. This test is modeled in the 2D axisymmetric case to achieve acceptable results. Application of a 200 kPa confining pressure is simulated as a uniformly distributed load applied normally inward to all external elements around the specimen; then, the cyclic load (q = 60 kPa & 70 kPa) is applied to the top of the sample.
In the Mohr-Coulomb model, since the test is carried out in an undrained condition (i.e., an undrained cyclic triaxial test), the internal friction angle of soil ( ) is zero and the cohesion value ( ) equals the undrained shear strength of soil ( ). Figure 7 shows the comparison of deviatoric stress versus axial strain in numerical simulations (kinematic hardening model and Mohr-Coulomb) with the experimental results in stress control mode. As shown in the diagrams, the behavior of undrained soils using the kinematic hardening model is in relatively good agreement with the results obtained from the experimental test, while the Mohr-Coulomb model fails to create the hysteresis loops, and the results thus differ from the experimental results. In addition, the material behavior in the Mohr-Coulomb model is in the elastic range. It is worth noting that this numerical simulation for the Mohr-Coulomb model is performed in two modes without considering and taking into account the damping ratios, and the results of the two cases are similar, which indicates that the linear elastic-perfectly plastic Mohr-Coulomb model cannot accurately predict the soil behavior under cyclic loading. In contrast, the kinematic model is capable of creating hysteresis loops and considers the soil stiffness degradation as the number of loading cycles increases. As shown in Figure 7, the undrained Young's modulus ratio of the numerical model to the experimental model ( ) increases with an increasing number of loading cycles, so that in the last cycle, it reaches its maximum value of approximately 1.4. Thus, a numerical model generally overestimates the soil stiffness degradation. The experimental curves in Figure 7 are related to tests C7 and C8.
Finally, after calibrating and using the corresponding equations, the values of the different hardening parameters are calculated for all five soil layers (see Figure 8). Thus, the values of the constitutive model parameters used in all of the analyses can be calibrated according to the second method (Table 1). It is worth noting that this numerical simulation for the Mohr-Coulomb model is performed in two modes without considering and taking into account the damping ratios, and the results of the two cases are similar, which indicates that the linear elastic-perfectly plastic Mohr-Coulomb model cannot accurately predict the soil behavior under cyclic loading. In contrast, the kinematic model is capable of creating hysteresis loops and considers the soil stiffness degradation as the number of loading cycles increases. As shown in Figure 7, the undrained Young's modulus ratio of the numerical model to the experimental model ( ) increases with an increasing number of loading cycles, so that in the last cycle, it reaches its maximum value of approximately 1.4. Thus, a numerical model generally overestimates the soil stiffness degradation. The experimental curves in Figure 7 are related to tests C7 and C8.
Finally, after calibrating and using the corresponding equations, the values of the different hardening parameters are calculated for all five soil layers (see Figure 8). Thus, the values of the constitutive model parameters used in all of the analyses can be calibrated according to the second method (Table 1).

Model Descriptions
As previously mentioned, many factors-such as the type and thickness of soil, shape, dimensions, and depth of a tunnel, as well as earthquake intensity-affect the seismic response of a tunnel to an earthquake. In this study, we aim to examine the effect of these factors on the seismic behavior of a tunnel and compare the results with those from previous studies performed using a linear elastic-perfectly plastic Mohr-Coulomb model [29]. In all of the studied models, five layers of soft soil are used with different characteristics and various thicknesses of 29.2 m, 19.6 m, 16 m, 7.2 m, and 3 m, corresponding to structural layers from bedrock to the soil surface, respectively. The schematic configuration of the models, including the soil layers and tunnel used in all of the analyses, is given in Figure 8, and the properties of the clay soils used in this study are summarized in Table 2.

Model Descriptions
As previously mentioned, many factors-such as the type and thickness of soil, shape, dimensions, and depth of a tunnel, as well as earthquake intensity-affect the seismic response of a tunnel to an earthquake. In this study, we aim to examine the effect of these factors on the seismic behavior of a tunnel and compare the results with those from previous studies performed using a linear elastic-perfectly plastic Mohr-Coulomb model [29]. In all of the studied models, five layers of soft soil are used with different characteristics and various thicknesses of 29.2 m, 19.6 m, 16 m, 7.2 m, and 3 m, corresponding to structural layers from bedrock to the soil surface, respectively. The schematic configuration of the models, including the soil layers and tunnel used in all of the analyses, is given in Figure 8, and the properties of the clay soils used in this study are summarized in Table 2. A circular tunnel with a diameter of 6 m and a square tunnel with a length of 5.3 m are used to investigate the effect of the shape of the tunnel on the dynamic response. Different tunnel-lining thicknesses of the circular tunnel are also considered to study the effect of lining thickness on the seismic response. The properties of the concrete used in the tunnel lining are given in Table 3. One of the factors affecting the seismic behavior of a tunnel is the intensity of the earthquake force imposed on the soil-tunnel model. The earthquake load applied to the base of models is that of the 1995 Kobe earthquake, i.e., 0.83 g PGA (peak ground acceleration) (as recorded at KJMA station); then, this is scaled down to 0.4 g to investigate the effect of this parameter on the dynamic soil-tunnel behavior. There are various methods [35][36][37][38][39] used for the boundary conditions of models to ensure that these conditions prevent the return of earthquake waves, and we use the viscous spring boundary method to surround our soil models. The acceleration time history of the Kobe earthquake, and the power amplitude and Fourier amplitude spectra used in this study, are presented in Figure 9. A circular tunnel with a diameter of 6 m and a square tunnel with a length of 5.3 m are used to investigate the effect of the shape of the tunnel on the dynamic response. Different tunnel-lining thicknesses of the circular tunnel are also considered to study the effect of lining thickness on the seismic response. The properties of the concrete used in the tunnel lining are given in Table 3. One of the factors affecting the seismic behavior of a tunnel is the intensity of the earthquake force imposed on the soil-tunnel model. The earthquake load applied to the base of models is that of the 1995 Kobe earthquake, i.e., 0.83 g PGA (peak ground acceleration) (as recorded at KJMA station); then, this is scaled down to 0.4 g to investigate the effect of this parameter on the dynamic soil-tunnel behavior. There are various methods [35][36][37][38][39] used for the boundary conditions of models to ensure that these conditions prevent the return of earthquake waves, and we use the viscous spring boundary method to surround our soil models. The acceleration time history of the Kobe earthquake, and the power amplitude and Fourier amplitude spectra used in this study, are presented in Figure  9.

Results and Discussions
In this study, the linearly elastic-perfectly plastic model and the kinematic hardening constitutive model are used to model soil behavior under earthquake load and some of the parameters that affect the seismic behavior of tunnels are investigated. The results that we obtain

Results and Discussions
In this study, the linearly elastic-perfectly plastic model and the kinematic hardening constitutive model are used to model soil behavior under earthquake load and some of the parameters that affect the seismic behavior of tunnels are investigated. The results that we obtain from a nonlinear kinematic hardening model are compared with those of a previous study (performed using a Mohr-Coulomb model) [29]. It is noteworthy that in this research, all the soil-tunnel models were simulated in addition to the kinematic hardening constitutive model using the Mohr-Coulomb model, and the results were in good agreement with the results of a previous study that considered the Mohr-Coulomb model using PLAXIS software. The results we obtain are discussed below.
To investigate the effect of tunnel embedment depth on the seismic response of a soil-tunnel system, a circular tunnel with a diameter of 6 m (D = 6 m) at depths of 6, 12, and 18 m (C1 = 6 m, C2 = 12 m, and C3 = 18 m) is subject to Kobe earthquake force with a maximum acceleration of 0.4 g to examine the effect of the interaction between soil and tunnel. Figure 10 shows the maximum acceleration response of the earthquake to the soil depth; as can be seen, the earthquake acceleration response decreases from the bedrock to the soil surface. Therefore, the forces that are exerted on deeper tunnels are larger. from a nonlinear kinematic hardening model are compared with those of a previous study (performed using a Mohr-Coulomb model) [29]. It is noteworthy that in this research, all the soiltunnel models were simulated in addition to the kinematic hardening constitutive model using the Mohr-Coulomb model, and the results were in good agreement with the results of a previous study that considered the Mohr-Coulomb model using PLAXIS software. The results we obtain are discussed below.
To investigate the effect of tunnel embedment depth on the seismic response of a soil-tunnel system, a circular tunnel with a diameter of 6 m (D = 6 m) at depths of 6, 12, and 18 m (C1 = 6 m, C2 = 12 m, and C3 = 18 m) is subject to Kobe earthquake force with a maximum acceleration of 0.4 g to examine the effect of the interaction between soil and tunnel. Figure 10 shows the maximum acceleration response of the earthquake to the soil depth; as can be seen, the earthquake acceleration response decreases from the bedrock to the soil surface. Therefore, the forces that are exerted on deeper tunnels are larger. The bending moment created in the tunnel lining is also investigated to determine the effect of tunnel depth on the seismic response. Figure 11 shows the bending moment created in a node on the lining at three different depths (6 m, 12 m, and 18 m) by considering the kinematic hardening model. The bending moment increases with increasing tunnel depth, which indicates an increase in the forces applied to the tunnel lining, i.e., the earthquake force increases with depth (see Figure 10). Figure 12 shows the comparison of tunnel displacement (considering the top node of the tunnel lining) at three different depths during the earthquake loading. As the tunnel depth increases, its displacement decreases, which means that the shallow tunnel experiences the largest displacement. The bending moment created in the tunnel lining is also investigated to determine the effect of tunnel depth on the seismic response. Figure 11 shows the bending moment created in a node on the lining at three different depths (6 m, 12 m, and 18 m) by considering the kinematic hardening model. The bending moment increases with increasing tunnel depth, which indicates an increase in the forces applied to the tunnel lining, i.e., the earthquake force increases with depth (see Figure 10). Figure 12 shows the comparison of tunnel displacement (considering the top node of the tunnel lining) at three different depths during the earthquake loading. As the tunnel depth increases, its displacement decreases, which means that the shallow tunnel experiences the largest displacement.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 23 Figure 11. Bending moment created at a node on the tunnel lining. Figure 11. Bending moment created at a node on the tunnel lining.
The acceleration response curve on the soil surface (at the center line of the model) for a tunnel at different depths (C1 = 6 m, C2 = 12 m, and C3 = 18 m) is shown in Figure 13, illustrating the effect of tunnel depth on the acceleration response of the soil surface. This graph results from a model in which the thickness of the tunnel lining is 0.3 m and the input motion of earthquake is 0.4 g. According to Figure 13, the maximum acceleration response on the soil surface does not change with increasing depth, and the value of maximum acceleration for all depths is approximately 0.06 g. Figure 11. Bending moment created at a node on the tunnel lining. The acceleration response curve on the soil surface (at the center line of the model) for a tunnel at different depths (C1 = 6 m, C2 = 12 m, and C3 = 18 m) is shown in Figure 13, illustrating the effect of tunnel depth on the acceleration response of the soil surface. This graph results from a model in which the thickness of the tunnel lining is 0.3 m and the input motion of earthquake is 0.4 g. According to Figure 13, the maximum acceleration response on the soil surface does not change with increasing depth, and the value of maximum acceleration for all depths is approximately 0.06 g.   The acceleration response curve on the soil surface (at the center line of the model) for a tunnel at different depths (C1 = 6 m, C2 = 12 m, and C3 = 18 m) is shown in Figure 13, illustrating the effect of tunnel depth on the acceleration response of the soil surface. This graph results from a model in which the thickness of the tunnel lining is 0.3 m and the input motion of earthquake is 0.4 g. According to Figure 13, the maximum acceleration response on the soil surface does not change with increasing depth, and the value of maximum acceleration for all depths is approximately 0.06 g. Figure 13. Acceleration response at the soil surface for a tunnel at various embedment depths, according to the kinematic hardening model. Three lining thicknesses (t1 = 0.15 m, t2 = 0.30 m, and t3 = 0.45 m) are considered for a circular tunnel of diameter 6 m to study the effect of the tunnel-lining thickness on its seismic behavior under the Kobe earthquake load. Figure 14 is based on the bending moment (normalized to the radius of the tunnel and the free-field shear stress ( M τR 2 ), (τ = Gγ max )) in terms of the lining thickness ratio ( t D ), and it shows a comparison of the results of this study using the nonlinear kinematic hardening and Mohr-Coulomb models with those of a previous study [29]. Based on these results, it can be seen that the tunnel-lining thickness affects the bending moment that occurs in the lining because the relative stiffness (i.e., the lining stiffness relative to the surrounding soil stiffness) varies with changes in the lining thickness. This means that the bending moment increases with the increasing lining thickness, indicating that the lining is able to withstand greater forces when it has a greater thickness. The earthquake forces of 0.4 g and 0.83 g are applied to the soil-tunnel system with a tunnel depth and thickness of 12 m and 0.3 m, respectively, to investigate the effect of the earthquake intensity on soil behavior. The distribution of equivalent plastic strain created near the tunnel for both the Mohr-Coulomb and kinematic hardening models is shown in Figure 15. Under the same conditions, the maximum equivalent plastic strain created in the soil elements in contact with the lining in the Mohr-Coulomb model is more than the kinematic hardening model, because the Mohr-Coulomb model cannot incorporate the damping caused by the soil plasticity. This means that the maximum value in the kinematic hardening model under 0.83 g (ε p = 0.0143) is approximately equal to the Mohr-Coulomb model under 0.4 g (ε p = 0.0126), but in the Mohr-Coulomb model, the equivalent plastic strain is zero in the soil elements away from the tunnel lining. Therefore, the Mohr-Coulomb model is not suitable for modeling soil under seismic loading, as it fails to generate loading and unloading loops, and it is unable to show the actual behavior of the soil. For instance, the equivalent plastic strain created during the earthquake loading (0.4 g and 0.83 g) for a node on the soil (above the tunnel crown) for both the Mohr-Coulomb and kinematic models is shown in Figure 16.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 23 Figure 13. Acceleration response at the soil surface for a tunnel at various embedment depths, according to the kinematic hardening model.
Three lining thicknesses (t1 = 0.15 m, t2 = 0.30 m, and t3 = 0.45 m) are considered for a circular tunnel of diameter 6 m to study the effect of the tunnel-lining thickness on its seismic behavior under the Kobe earthquake load. Figure 14 is based on the bending moment (normalized to the radius of the tunnel and the free-field shear stress ( 2 ), ( = )) in terms of the lining thickness ratio ( t ⁄ ), and it shows a comparison of the results of this study using the nonlinear kinematic hardening and Mohr-Coulomb models with those of a previous study [29]. Based on these results, it can be seen that the tunnel-lining thickness affects the bending moment that occurs in the lining because the relative stiffness (i.e., the lining stiffness relative to the surrounding soil stiffness) varies with changes in the lining thickness. This means that the bending moment increases with the increasing lining thickness, indicating that the lining is able to withstand greater forces when it has a greater thickness.  , ⁄ = . and ⁄ = . .
The earthquake forces of 0.4 g and 0.83 g are applied to the soil-tunnel system with a tunnel depth and thickness of 12 m and 0.3 m, respectively, to investigate the effect of the earthquake intensity on soil behavior. The distribution of equivalent plastic strain created near the tunnel for both the Mohr-Coulomb and kinematic hardening models is shown in Figure 15. Under the same conditions, the maximum equivalent plastic strain created in the soil elements in contact with the lining in the Mohr-Coulomb model is more than the kinematic hardening model, because the Mohr-Coulomb model cannot incorporate the damping caused by the soil plasticity. This means that the maximum value in the kinematic hardening model under 0.83 g ( = 0.0143) is approximately equal to the Mohr-Coulomb model under 0.4 g ( = 0.0126), but in the Mohr-Coulomb model, the equivalent plastic strain is zero in the soil elements away from the tunnel lining. Therefore, the Mohr-Coulomb model is not suitable for modeling soil under seismic loading, as it fails to generate loading and unloading loops, and it is unable to show the actual behavior of the soil. For instance, the equivalent plastic strain created during the earthquake loading (0.4 g and 0.83 g) for a node on the soil (above the tunnel crown) for both the Mohr-Coulomb and kinematic models is shown in Figure  16.    The choice of the shape and dimensions of a tunnel depends on its intended purpose. For example, circular, quadrilateral and horseshoe cross-section shapes are used for railways tunnels, and each of these tunnels exhibits different behavior under seismic loads. Here, we study two circular and square tunnels with the same embedment depth (D = 12 m) and lining thickness (t = 0.30 m), which are subject to the same earthquake intensity (0.4 g), to investigate the effect of tunnel shape on the seismic behavior. The difference between the results we obtain and those from the previous study [29] (Figure 17) are due to the different constitutive model used. The choice of the shape and dimensions of a tunnel depends on its intended purpose. For example, circular, quadrilateral and horseshoe cross-section shapes are used for railways tunnels, and each of these tunnels exhibits different behavior under seismic loads. Here, we study two circular and square tunnels with the same embedment depth (D = 12 m) and lining thickness (t = 0.30 m), which are subject to the same earthquake intensity (0.4 g), to investigate the effect of tunnel shape on the seismic behavior. The difference between the results we obtain and those from the previous study [29] (Figure 17) are due to the different constitutive model used. The maximum bending moment generated at all nodes on the tunnel lining for both Mohr-Coulomb and kinematic hardening models are shown in Figure 17a,b, which are for the circular tunnel and square tunnel, respectively. The horizontal axis of these graphs indicates the perimeter of the tunnel (zero corresponds to the top point of the tunnel, and other points are selected based on their clockwise-perimeter distance from the top point). Figure 18 shows the comparison of the bending moment created for one node on the top of the circular tunnel (tunnel crown) during the The maximum bending moment generated at all nodes on the tunnel lining for both Mohr-Coulomb and kinematic hardening models are shown in Figure 17a,b, which are for the circular tunnel and square tunnel, respectively. The horizontal axis of these graphs indicates the perimeter of the tunnel (zero corresponds to the top point of the tunnel, and other points are selected based on their clockwise-perimeter distance from the top point). Figure 18 shows the comparison of the bending moment created for one node on the top of the circular tunnel (tunnel crown) during the earthquake loading for both constitutive models. According to Figure 17a,b and Figure 18, there is a relatively large difference between the results of the two models, so that in the Mohr-Coulomb model, the absolute of the bending moment values at all nodes is more than that of the nonlinear kinematic hardening model. The results obtained from this study of the effect of the tunnel shape on the seismic behavior show that the maximum bending moment generated on the perimeter of the square tunnels is more than that generated on the perimeter of circular tunnel. It also shows that the corners of the square tunnel have maximum values, proving that the circular tunnel performs better than the square tunnel under the applied seismic loads.

Conclusions
The use of advanced constitutive models is limited in practice due to the difficulty in calibrating several constitutive parameters. Here, a linear elastic model is used for tunnels, and a linear elasticperfectly plastic Mohr-Coulomb model and a nonlinear kinematic hardening model are used as advanced constitutive models for clay soil. Two different methods-one using data from a cyclic shear test simulation and the experimental curve of Ishibashi and Zhang [32] and the other using experimental data and some data from a simulation of an undrained cyclic triaxial test (which is a new method)-are used to calibrate the parameters of a kinematic hardening constitutive model. The results of this research can be summarized as follows.
1. Choosing the appropriate constitutive model to achieve logical results is key to evaluating the seismic behavior of the tunnel in numerical methods. The use of a nonlinear kinematic hardening model that considers the effect of soil stiffness degradation is appropriate for structures under cyclic loads. 2. Calibrating the parameters of the kinematic hardening constitutive model using the data from undrained cyclic triaxial tests afford results that are almost identical to those from using the Ishibashi and Zhang equations and data from simulation of the cyclic shear tests. 3. The nonlinear kinematic hardening constitutive model can better demonstrate the cyclic deformation behavior of soils under cyclic loading than conventional and simple models such as the linear elastic-perfectly plastic Mohr-Coulomb model. 4. The Mohr-Coulomb model overestimates the design and it is not economically acceptable. For example, the bending moments created in the tunnel lining are much larger in the Mohr-Coulomb model than in the kinematic model. 5. The plastic strain of soil increases in both the kinematic hardening and Mohr-Coulomb models as the intensity of the earthquake increases from 0.4 to 0.83 g, but this increase is greater in the Mohr-Coulomb model due to its inability to create of hysteresis loops. Therefore, the Mohr-

Conclusions
The use of advanced constitutive models is limited in practice due to the difficulty in calibrating several constitutive parameters. Here, a linear elastic model is used for tunnels, and a linear elastic-perfectly plastic Mohr-Coulomb model and a nonlinear kinematic hardening model are used as advanced constitutive models for clay soil. Two different methods-one using data from a cyclic shear test simulation and the experimental curve of Ishibashi and Zhang [32] and the other using experimental data and some data from a simulation of an undrained cyclic triaxial test (which is a new method)-are used to calibrate the parameters of a kinematic hardening constitutive model. The results of this research can be summarized as follows.

1.
Choosing the appropriate constitutive model to achieve logical results is key to evaluating the seismic behavior of the tunnel in numerical methods. The use of a nonlinear kinematic hardening model that considers the effect of soil stiffness degradation is appropriate for structures under cyclic loads.

2.
Calibrating the parameters of the kinematic hardening constitutive model using the data from undrained cyclic triaxial tests afford results that are almost identical to those from using the Ishibashi and Zhang equations and data from simulation of the cyclic shear tests.

3.
The nonlinear kinematic hardening constitutive model can better demonstrate the cyclic deformation behavior of soils under cyclic loading than conventional and simple models such as the linear elastic-perfectly plastic Mohr-Coulomb model.

4.
The Mohr-Coulomb model overestimates the design and it is not economically acceptable. For example, the bending moments created in the tunnel lining are much larger in the Mohr-Coulomb model than in the kinematic model.

5.
The plastic strain of soil increases in both the kinematic hardening and Mohr-Coulomb models as the intensity of the earthquake increases from 0.4 to 0.83 g, but this increase is greater in the Mohr-Coulomb model due to its inability to create of hysteresis loops. Therefore, the Mohr-Coulomb model can't be used to accurately predict the behavior of soil under earthquake loading. 6.
Changing the depth of a tunnel has no effect on the maximum acceleration response on the soil surface, but as the depth increases, greater forces are applied to the lining, and the bending moment created on the tunnel lining also increases. The lining displacement also decreases with the increasing tunnel depth. 7.
Another parameter that affects the dynamic behavior of tunnels is the tunnel-lining thickness. The bending moment in the tunnel lining increases with its thickness so that as the dimensions of the tunnel increase, the EI increases and the structure becomes more stiff and able to withstand more forces. The resistance to deformation and the flexibility also increase under dynamic loads. 8.
The shape of the tunnel dictates the overall response of the soil-tunnel system under dynamic loads. Circular tunnels show better performance than square tunnels against seismic loads. Generally, rounding the corners in square tunnels causes the tunnels to perform better under the imposed loads.
Author Contributions: M.S.A. and X.C. equally contributed to the conceptualization, analysis, and methodology. M.S.A. wrote the original draft; X.C. contributed mainly to the conceptualization, review, and editing. All authors have read and agreed to the published version of the manuscript.
Funding: The first author of this article is financially supported by Chines Government Scholarship (CSC).

Conflicts of Interest:
The authors declare no conflicts of interest.