A stable CS-FEM for the static and seismic stability of a single square tunnel in the soil where the shear strength increases linearly with depth

H.C. Nguyen, L. Nguyen-Son

a Department of Civil and Environmental Engineering, Imperial College London, London, UK

b Department of Civil Engineering and Industrial Design, University of Liverpool, Liverpool, UK

c Faculty of Civil Engineering, Ho Chi Minh City University of Technology (HUTECH), Ho Chi Minh City, Viet Nam

Keywords:Tunnels Stability Limit analysis Cell-based smoothed finite element method(CS-FEM)Second-order cone programming (SOCP)

ABSTRACT A numerical procedure using a stable cell-based smoothed finite element method(CS-FEM)is presented for estimation of stability of a square tunnel in the soil where the shear strength increases linearly with depth. The kinematically admissible displacement fields are approximated by uniform quadrilateral elements in conjunction with the strain smoothing technique, eliminating volumetric locking issues and the singularity associated with the Mohr-Coulomb model.First,a rich set of simulations was performed to compute the static stability of a square tunnel with different geometries and soil conditions. The presented results are in excellent agreement with the upper and lower bound solutions using the standard finite element method (FEM). The stability charts and tables are given for practical use in the tunnel design,along with a newly proposed formulation for predicting the undrained stability of a single square tunnel. Second, the seismic stability number was computed using the present numerical approach.Numerical results reveal that the seismic stability number reduces with an increasing value of the horizontal seismic acceleration (αh), for both cases of the weightless soil and the soil with unit weight. Third, the link between the static and seismic stability numbers is described using corrective factors that represent reductions in the tunnel stability due to seismic loadings. It is shown from the numerical results that the corrective factor becomes larger as the unit weight of soil mass increases;however, the degree of the reduction in seismic stability number tends to reduce for the case of the homogeneous soil. Furthermore, this advanced numerical procedure is straightforward to extend to three-dimensional (3D) limit analysis and is readily applicable for the calculation of the stability of tunnels in highly anisotropic and heterogeneous soils which are often encountered in practice.

1. Introduction

In this paper, we present a cell-based smoothed finite element method (CS-FEM) to estimate the static and seismic stability of a square tunnel in the soil where the shear strength varies with depth in a linearly increasing manner. Sloan and Assadi (1991)presented two numerical procedures using finite element formulations of the classic lower and upper bound theorems to calculate the rigorous bounds on the undrained stability of a square tunnel.The loads to resist the collapse(Sloan and Assadi,1991)were solved by optimizations established in large sparse linear programming problems (LNP) using limit theorems. An extension to the establishment of optimizations of two bounds in a conic form was outlined by Wilson et al. (2013), along with the use of rigid-block upper bound method to predict the undrained stability of the square tunnel. The solutions to the undrained stability of a square tunnel were improved significantly, with upper and lower bound solutions to the collapse load in a range of 5%. Recently, Vo-Minh and Nguyen-Son (2021) presented an upper bound limit analysis using the node-based smoothed finite element to calculate the stability of two circular tunnels at different depths in cohesivefrictional soils, with very satisfactory results. It was reported by Meng et al. (2020) and Nguyen and Vo-Minh (2022a) that such a low-order mixed element can overcome the volumetric locking problems as well as the singularity associated with the Mohr-Coulomb model. More recently, Nguyen and Vo-Minh (2022b)confirmed that the use of the strain smoothing technique in the limit analysis offers some advantages: (i) reductions in the size of optimization problem; and (ii) productions of stable and accurate plasticity solution to seismic bearing capacity. In this study, we exploit the advantages of the CS-FEM (Nguyen, 2021; Nguyen and Vo-Minh, 2022b) based on the quadrilateral element mesh to revisit the upper bound solution to the undrained stability of a square tunnel. Furthermore, several analyses were conducted to calculate the seismic stability of a single square tunnel by including the horizontal seismic acceleration αhin the simulations.This study is the extended work of Nguyen(2021)who computed the seismic stability of tunnel in homogeneous soil considering the case that the shear strength varies linearly with depth.

First, the single square tunnel problem was discretized using uniform rectangular element meshes, followed by a strain smoothing technique to approximate the displacement fields.Second, the upper bound limit analysis was cast as the secondorder cone programming (SOCP) problem. This numerical procedure was adopted to assess the stability of a square tunnel that is of a width B and rests at depth H, as illustrated in Fig.1. The collapse associated with the surcharge σsand the soil weight γ was used to determine the internal tunnel pressure σtneeded to resist the collapse. The tunnel is constructed in undrained conditions where the undrained shear strength linearly increases with depth:

where cu0and cu(z)are the cohesions at the ground and the depth z, respectively.

According to Sloan and Assadi (1991), the overall stability of tunnel can be represented by two load terms,i.e.(σs-σt)/cu0and γB/cu0,which significantly depend on the magnitudes of H/B and ρB/cu0.Wilson et al. (2013) focused on the so-called stability number as below:

Intensive investigations into how the undrained stability number Nsvaries with the depth and the rate of changes in the undrained shear strength with depth were performed by Wilson et al.(2013)using the combination of FEM and SOCP,producing a better set of accurate collapse loads when compared with the use of FEM and LNP in the upper bound procedure (Sloan and Assadi, 1991).This research extends the work of Sloan and Assadi (1991) and Wilson et al.(2013)by giving a new(better)upper bound solution to the undrained stability number Nsusing the numerical approach based on the CS-FEM and SOCP. Although the present numerical procedure is applied to resolving the stability of tunnels in purely cohesive soil,this upper bound analysis using CS-FEM is readily to estimate stability numbers for the case of highly anisotropic and heterogeneous soils which are widely encountered in the tunnel design by engineers (Krabbenhoft and Lyamin, 2015; Krabbenhoft et al., 2019). In addition, it is worth noting that the cohesion reduces as the displacement increases beyond the initial failure stage as noted by Zhang et al. (2017). This fact does not affect the magnitude of the stability number calculated using the limit theorem; however, the cohesion softening is of significance in analyzing the progressive failure of geotechnical problems,such as landslides (Skempton and Hutchinson,1969).

Fig. 1. Illustration of a square tunnel in the soil where the shear strength increases linearly with depth (under the static (αh = 0) and seismic conditions (αh ≠0)).

The paper is structured as follows. In Section 2.1, we shortly present how to apply the strain smoothing technique over a single smoothing cell using uniform quadrilateral mesh. The formulation of the upper bound procedure using CS-FEM and SOCP is subsequently given in Section 2.2 to compute both the static and seismic stability of a square tunnel in heterologous soils. Several simulations using different numerical methods (both the smoothed and the standard FEMs) were conducted to prove the efficiency and accuracy of the analyses.In addition,a rich set of simulations were intensively performed in Section 3; and the numerical results are then checked against the prior contributions (Sloan and Assadi,1991; Wilson et al., 2013). Section 4 includes some highlighting features that emerge from the present results.

2. Upper bound analysis using CS-FEM and SOCP

2.1. An overview of the CS-FEM

Fig. 2. Division of a quadrilateral element into smaller smoothing cells (SSC): (a) SSC = 1; (b) SSC = 2; (c) SSC = 3; (d) SSC = 4; (e) SSC = 8; and (f) SSC = 16.

where Γcis the boundary of smoothing domain as illustrated in Fig. 3, and the normal vector nxis expressed as

Therefore,the smoothing strain rate can be calculated as～.εxc=.B.d, in which d is the displacement vector of the nodes associated with the quadrilateral element, and .B is the strain-displacement matrix defined by

Fig. 3. The smoothing domain using a single smoothing cell Ωk with the Gauss point xkG of boundary segment Γkc.

2.2. Establishment of upper bound optimization problems as SOCP

It is well known that the structure will collapse if and only if there exists a kinematically admissible displacement field .u?U,where U is a space of kinematically admissible velocity field, thus we have

The space of kinematically admissible velocity field U is denoted by

For plane strain problems based on the Mohr-Coulomb failure criterion, Makrodimopoulos and Martin (2007) proposed the internal plastic dissipation equation as follows:

where c and φ are the cohesion and friction angle of the soil,respectively.

For an associated flow rule,the plastic strain rate vector is given by

where μ is a non-negative plastic multiplier, and the Mohr-Coulomb yield function ψ(σ) can be expressed in the form of stress components as

We aim to form the upper bound solutions as conic quadratic optimization in the quadratic form as noted by Zhang et al.(2019):

Using CS-FEM, the problem is discretized by Ncell smoothing domains.The smoothed strain rate～.ε can be calculated from Eq.(3).The upper bound limit analysis for plane strain problems using the Mohr-Coulomb failure criterion can be written as

subject to

where Aiis the area of the smoothing domain of cell i, and the collapse load multiplier β+is the static or seismic bearing capacity factor in this study. The last constraint in Eq. (19b) is expressed in the conic form.As a result,the conic interior-point optimizer of the academic MOSEK package(MOSEK ApS,2015)is used to solve this problem. The upper bound using the CS-FEM has been written using the Matlab language. The computations were performed on the Dell precision 5520 (Intel(R) Xeon(R) CPU E3-1505 M v5 2.80 GHz,32 Gb RAM)in the Window 10 Pro;and the combination of CS-FEM and SOCP in the upper bound analysis has been written in the Matlab language.In short,the whole procedure of numerical procedure using CS-FEM and SOCP can be depicted in Fig. 4 that shows algorithms to implement into the code.

3. Limit analysis of stability of a square tunnel

Having presented a combination of CS-FEM and SOCP in the previous section,we used this numerical procedure to calculate the so-called Ns= (σs-σt)/cu0for various H/B ratios. As shown in Fig. 5, only half of a square tunnel was considered in the analyses.The boundary conditions were set in the simulations as given in Fig. 5 due to the geometry of symmetry. Displacement conditions along the boundary were set in the simulations as illustrated in Fig. 5, where the displacement condition u = v = 0 was installed along the bottom and the right of domain considered and the horizontal displacements u=0 was set along the left of the domain of analyses. It is noted that all domains should be large enough to remove the effect of boundary conditions on the numerical solutions.

Fig. 4. Illustration of upper bound limit analysis using CS-FEM and SOCP.

Fig. 5. Representative quadrilateral element mesh for H/B = 3, showing boundary conditions.

Table 1 The static stability of a square tunnel for the case of H/B = 1, ρB/cu0 = 0, and γB/cu0 = 0.

In order to prove the accuracy and efficiency of the current numerical approach using CS-FEM,we consider the case of H/B=1,ρB/cu0= 0, and γB/cu0= 0. The total number of variables Nvar, the number of elements Ne,degree of freedoms,Sdof,CPU time,and the values of the static stability number are listed in Table 1.The lower and upper bound solutions to the stability number given by Wilson et al. (2013) are Ns= 1.94 and Ns= 1.98, respectively. While this study used only about more than 7000 elements, the present results are equal to the finite element upper bound limit analysis given by Wilson et al. (2013) who used a significantly larger number of elements in the analyses (i.e. about 100,000 elements with discontinuity elements). The validation of this numerical approach for the collapse loads under seismic conditions has been presented by Nguyen and Vo-Minh (2022b), and the interested readers are referred to this reference for more details on the comparison of the use of CS-FEM and other FEMs in the limit analysis. This confirms the effectiveness and the accuracy of the present analysis for calculation of the stability of a square tunnel.

It is worth noting that the collapse of the tunnel occurs at constant volume with presented simulations using the constrained conditions as stated by Wilson et al. (2013) (see Eq. (19b)). The dimensionless stability is quantified by setting the surcharge to zero, i.e. σs= 0. A rich set of simulations were performed by the present approach for various cases of dimensionless parameters(H/B,ρB/cu0,and γB/cu0).We varied H/B from 1 to 10 to investigate how the depth of a square tunnel influences the stability.For each tunnel depth,the soil unit weight(γB/cu0)was varied from 0 to 5,in which γB/cu0= 0 is analogous to the case of tunnel in the weightless soil.In addition, the undrained shear strength varies linearly with depth, representing the strength parameter (ρB/cu0) of the soil changing between 0 and 1. The dimensionless parameter Nsis computed and reported in the Appendix,along with a comparison of the present results with two bounds on the undrained stability given by Sloan and Assadi(1991)and Wilson et al.(2013),as shown in Figs. 6-10. Overall, the numerical results are better than both upper bounds given in Wilson et al. (2013) and the upper bound solutions reported in Sloan and Assadi (1991). For all cases of analyses,the stability numbers obtained using CS-FEM are lower than those given by the prior contributions (Sloan and Assadi, 1991;Wilson et al., 2013). As expected, as the shear strength increases,the stability number increases for all cases of analysis of the weightless soil (γB/cu0= 0) and the soil considering its weight(γB/cu0= 1-5).As noted by Wilson et al.(2013),a negative value of the stability number indicates that the application of a compressive normal stress to the tunnel face is required to prevent the failure from occurring.On the other hand,a positive value of Nsimplies that the tunnel itself is safe from failure.Numerical results reveal that, for a tunnel of moderate depth, Nsvalues decrease significantly as γB/cu0increases, meaning that the supporting forces are needed to maintain the stability of the tunnel.

Fig. 6. Comparison of the stability of a square tunnel obtained in the present study and those given by Wilson et al. (2013) for: (a) H/B = 1; and (b) H/B = 2.

Fig. 7. Comparison of the stability of a square tunnel obtained in the present study and those given by Wilson et al. (2013) for: (a) H/B = 3; and (b) H/B = 4.

Fig. 8. Comparison of the stability of a square tunnel obtained in the present study and those given by Wilson et al. (2013) for: (a) H/B = 5; and (b) H/B = 6.

Fig. 9. Comparison of the stability of a square tunnel obtained in the present study and those given by Wilson et al. (2013) for: (a) H/B = 7; and (b) H/B = 8.

Fig.10. Comparison of the stability of a square tunnel obtained in the present study and those given by Wilson et al. (2013) for: (a) H/B = 9; and (b) H/B = 10.

For practical design,three rigid block mechanisms(Wilson et al.,2013)are useful to calculate the stability number of a single square tunnel.It was reported by Wilson et al.(2013)that the mechanism a is suitable for estimations of the stability number for a shallow tunnel where the failure occurs on the top (i.e. the roof) of tunnel(Fig. 11). Numerical results are found in an excellent agreement with the trapdoor style mechanism (the mechanism a) for the shallow tunnel, as shown in Fig.11. When the ratio H/B increases,the use of the mechanisms b and c are likely to successfully capture the failure mechanism which takes place on the wall and the floor of tunnel. Figs. 12 and 13 show the numerical results of failure mechanism for deeper tunnels, agreeing very well with the mechanisms b and c proposed by Wilson et al.(2013).In particular,Fig. 12 shows the comparison of failure mechanism with the mechanism b for H/B = 1, ρB/cu0= 0.25, and γB/cu0= 3, in which both types of failure occur on the wall of the tunnel. For deeper tunnels, the failure domain becomes deeper, as shown in Fig. 13,where the failure occurs beneath the floor of tunnel. This observation provides a check on applications of the rigid block mechanism c for prediction of the undrained stability of a deeper tunnel.

Fig. 11. A comparison of the power dissipation intensity and the rigid block mechanism for H/B = 1, ρB/cu0 = 0, and γB/cu0 = 0: (a) This study; and (b) The mechanism a given in Wilson et al. (2013).

Fig.12. A comparison of the power dissipation intensity and the rigid block mechanism for H/B = 1, ρB/cu0 = 0.25, and γB/cu0 = 3: (a) This study; and (b) The mechanism b given in Wilson et al. (2013).

In addition, the magnitude of the stability number for the homogeneous soil obtained by the present analysis is better than both upper bound solutions given in Sloan and Assadi(1991)and Wilson et al.(2013).It was revealed by Wilson et al.(2013)that the failure pattern is strongly dependent upon the H/B and γB/cu0ratios. A comparison of power dissipation intensity with the work of Wilson et al.(2013)for H/B=1,ρB/cu0=0,and γB/cu0=0 is given in Fig.14,confirming that the trapdoor mechanism successfully captures the failure pattern for shallow tunnel in the soil with a low unit weight.Numerical results confirm that the failure pattern of a shallow tunnel with H/B ≤1 in a weightless soil remains unchanged,having a trapdoor mechanism when the ratio ρB/cu0varies.However,this feature does not hold for a shallow circle tunnel as reported by Wilson et al.(2011).In addition,a rigid block mechanism(b)seems to be suitable for capturing the failure mode for a shallow tunnel(H/B ≤1) in the soil considering its weight.

Fig.13. A comparison of the power dissipation intensity and the rigid block mechanism for H/B = 9, ρB/cu0 = 0, and γB/cu0 = 0: (a) This study; and (b) The mechanism c given inWilson et al. (2013).

Fig. 14. Comparison of power dissipation intensity with the work of Wilson et al.(2013)for H/B=1,ρB/cu0=0,and γB/cu0=0:(a)Wilson et al.(2013);and(b)CS-FEM.

Fig.15. Comparison of the power dissipation intensity for H/B=4,ρB/cu0= 0,and γB/cu0 = 1: (a) Wilson et al. (2013); and (b) CS-FEM.

Figs. 15-18 show the comparison of the failure mechanism obtained in this study with those given by Wilson et al. (2013)when the depth of tunnels increases, indicating that the failure shapes derived from both numerical approaches are in an excellent agreement. The expansion of the domain of failure mechanism in conjunction with the depth of failure zone provides a check on observations given in Wilson et al.(2013).In addition,as the ρB/cu0ratio increases, the failure becomes more localized when comparing the failure domain in Fig. 15 with that in Fig. 16 for a tunnel of moderate depth (H/B = 4). The effect of ρB/ cu0on characteristics of failure pattern becomes significant as the H/B ratio increases,as illustrated in Figs.17 and 18 for γB/cu0= 1.Although the failure pattern becomes narrower with increasing values of ρB/cu0, the value of stability number increases, as reported in the Appendix.

Fig.16. Comparison of the power dissipation intensity for H/B=4,ρB/cu0=1,and γB/cu0 = 1: (a) Wilson et al. (2013); and (b) CS-FEM.

To facilitate the practical use of tunnel design by engineers,we propose a new expression to compute the undrained stability of a square tunnel in terms of the three dimensionless parameters, i.e.H/B,ρB/cu0,and γB/cu0.It is shown in Figs.6-10 where the so-called Ns= (σs-σt)/cu0is proportional to the γB/cu0ratio in a linear manner,leading to a formulation of the stability number that is the function of N0and Nγ:

where N0and Nγ are the stability factors which account for the weightless soil and the soil considering weight, respectively. It is noted that two parameters, ζ and ψ, are quantified using the appropriate regression over the numerical results.Fitting the whole set of numerical results listed in the Appendix for the case of γB/cu0=0 gives a new expression for N0which is a function of H/B and ρB/cu0:

Fig.17. Comparison of the power dissipation intensity for H/B=7,ρB/cu0= 0,and γB/cu0 = 3: (a) Wilson et al. (2013); and (b) CS-FEM.

Fig.18. Comparison of the power dissipation intensity for H/B=7,ρB/cu0=1,and γB/cu0 = 3: (a) Wilson et al. (2013); and (b) CS-FEM.

In order to generate an expression for the factor Nγ, various options using parametric curve-fitting techniques to fit Eqs. (20)and (21) are adopted until an appropriate expression can be achieved to match quite closely with the present solution. The following is the parametric equation for the stability factor accounting for the soil weight:where the two parameters ζ and ψ(Eq.(20))are the functions of H/B.A check on the newly proposed equation was made by comparing the values of Nsderived from Eq.(20)and the equation proposed by Wilson et al. (2013), along with the average values of both lower and upper bounds on the Nsvalues. As shown in Fig.19, the new parametric equation gives much better values of the undrained stability of a circular tunnel than that formulation given by Wilson et al.(2013).For the weightless soil,the values of Nsare equivalent to the case of N0.It can be observed from Fig.19a that there are no differences between the values obtained by the new equation and the average values of two bounds given by Wilson et al.(2013)for the homogeneous soil.In particular,both equations give almost the same values of the stability number for ρB/cu0≤0.25.However,it is shown that the new equation derived from the present study better predicts the undrained stability for ρB/cu0≥0.5.Similar behavior is obtained for the cases of γB/cu0=1 for ρB/cu0≤0.25.On the other hand, when ρB/cu0increases, the present formulation predicts better values(i.e.closer to the average values of two bounds on Ns)of the undrained stability of a square tunnel when compared with that given by Wilson et al. (2013). For the soil considering weight,when γB/cu0varies from 1 to 5, the Nsvalues obtained using the new expression are much closer to the average values of the two bound solutions, indicating that the new formulation gives better values of the stability number than that given by Wilson et al.(2013). It is shown in Fig.19b-f that for deeper tunnels, the new closed-form expression of Nsis much better than that given by Wilson et al.(2013)when utilized to predict the stability number.In short,for both the weightless soil and the soil considering weight,the new equation derived from the present study acts as a better means to interpret the undrained stability of a square tunnel.Although the new equation proposed for predicting the undrained stability number gives a better fit to the exact values of the stability number for the purely cohesive soil, there are still trivial discrepancies between the predicted results from the current expression and the average results for highly anisotropic and heterogeneous soils (i.e. with large values of ρB/cu0). These predictions can be applicable in the practical use by engineers.

Apart from calculation of static stability of a square tunnel,several simulations were carried out to estimate the seismic stability of a square tunnel. Herein, an extension of the work carried out by Nguyen(2021),who assessed the effect of seismic loading on the tunnel stability for a homogeneous soil (i.e. ρB/cu0= 0), was made by considering the case that the undrained shear strength linearly increases with depth. It is shown from Fig. 20 that under the static condition,the failure mechanism is symmetric;however,the plastic mechanism becomes asymmetric with the increasing value of horizontal seismic acceleration.This leads to a reduction in the magnitude of the stability number.

In order to assess the effect of seismic loadings on the undrained stability of tunnel in the soil where the shear strain increases linearly with depth, several simulations were performed for two values of depth ratio,H/B=1 and H/B=3.We varied αhfrom 0 to 0.5,taking into account the weightless soil(γB/cu0=0)and the soil considering weight (γB/cu0=1) in the analyses for each depth. In this case, the so-called seismic stability number is defined as

4. Conclusions

A CS-FEM is presented to calculate the undrained stability of a square tunnel, eliminating volumetric locking issues and the singularity associated with the Mohr-Coulomb model. A series of computations was performed to calculate the undrained stability of a square tunnel in the soil where the undrained shear strength varies with depth in a linear manner. Numerical results are in excellent agreement with prior contributions (Sloan and Assadi,1991; Wilson et al., 2013), giving new upper bounds on the undrained stability that are better than the results given by Sloan and Assadi(1991)and Wilson et al.(2013)for deeper tunnels.As γB/cu0increases, the numerical results of Ns= (σs-σt)/cu0become significantly improved. Failure mechanisms obtained agree very well with the rigid block mechanisms used in Sloan and Assadi(1991) and the kinematic failure patterns used in the upper bound FEM, which changes with increasing depth values. Several stability tables and charts are given for the practical use of tunnel design, along with the newly proposed approximate equation for calculating the undrained stability of a square tunnel in the soil where the shear strength increase with depth in a linear form.

Fig.19. Comparison of the N0 and Nγ values obtained by the newly proposed expression with those given by Wilson et al.(2013):(a)N0 for γB/cu0=0;(b)Nγ for γB/cu0=1;(c)Nγ for γB/cu0 = 2; (d) Nγ for γB/cu0 = 3; (e) Nγ for γB/cu0 = 4; and (f) Nγ for γB/cu0 = 5. LB and UB denote the lower and upper bounds, respectively.

Several computational performances have been conducted to calculate the seismic stability number.The results confirm that the seismic stability number reduces as αhincreases for both cases of the weightless soil and the soil considering weight. The corrective coefficient was quantified to link the static stability number to its seismic counterpart. The numerical results reveal that the corrective coefficient esEbecomes larger with increasing values of soil unit weight. In addition, the magnitude of corrective coefficient increases as the γB/cu0ratio reduces; however, this coefficient drops significantly for deep tunnels (i.e. the H/B ratio increases).This upper bound numerical procedure is readily to extend to three-dimensional (3D) limit analysis of tunnel stability in highly anisotropic and heterogeneous soils which are often encountered in tunnel designs by geotechnical engineers.

Fig. 20. Displacement field and power dissipation intensity for the static and seismic stability of a square tunnel with γB/cu0 = 3, ρB/cu0 = 3, and H/B = 3: (a) power dissipation intensity for Ns=5.87 and αh=0;(b)Displacement field for Ns=5.87 and αh=0;(c)Power dissipation intensity for Ns=5.72 and αh=0.1;and(d)Displacement field for Ns=5.72 and αh = 0.

Table 2 The seismic stability number and the corrective coeffcient of a square tunnel for H/B = 1, and γB/cu0 = 0 and 1, with A=ρB/cu0.

Table 3 The seismic stability number and the corrective coeffcient of a square tunnel for H/B = 3, and γB/cu0 = 0 and 1, with A =ρB/cu0.

Table 4 The seismic stability number and the corrective coeffcient of a square tunnel for H/B = 3, and γB/cu0 = 2 and 3, with A =ρB/cu0.

Table 5 The seismic stability number and the corrective coeffcient of a square tunnel for H/B = 3, and γB/cu0 = 4 and 5, with A=ρB/cu0.

Fig.21. The seismic stability numberand the corrective coeffcient for H/B=1:(a)N for γB/cu0 = 1;(b) for γB/cu0 = 0;(c)h for γB/cu0 = 1;and(d) for γB/cu0 = 0.

Fig.22. The seismic stability number and the corrective coeffcient for H/B=3:(a)Nfor γB/cu0 = 1;(b) for γB/cu0 = 0;(c)e for γB/cu0 = 1;and(d)e for γB/cu0 = 0.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This is part of the TPS project. The first author is grateful to a Vied-Newton PhD scholarship and a Dixon scholarship from Imperial College London, UK, for supporting his studies at Imperial College London. He is also indebted to the Dean’s Fund from Imperial College London for financial support (2017-2020).