The stability evaluation index for large-section clay tunnels is quantitative, and it can estimate the variability between the mechanical state of the surrounding rock after tunnel excavation and its limit state.

(1) The stability reserve factor for the large cross-section clay tunnel

Firstly, the actual mechanical state of the surrounding rock at one point after tunnel excavation is S ₀ and the limit state at that point is S _L , by analogy with the slope safety factor, the calculated expression for the tunnel stability reserve factor F _s is given as follows. F s = S L S 0 (1)

The tunnel stability reserve factor can be applied to different geological conditions, environmental conditions, additional loading factors, et al., it only depends on the actual mechanical state and the limit state of the surrounding rock and does not depend on the constitutive model. Therefore, it is well suited for the stability evaluation of a large cross-section clay tunnel.

Figure 1 gives the correlation between the tunnel stability reserve factor based on Eq. 1 and the actual mechanical state of the surrounding rock. The ultimate state of the surrounding rock is closely related to the physico-mechanical properties of the clay. The process of deterioration of the physico-mechanical properties of the clay is specifically characterized as a change in the limit state of the surrounding rock. When the physico-mechanical properties of the surrounding rock remain unchanged, then the actual mechanical state of the surrounding rock can be adjusted to meet a given tunnel stability reserve factor.

(2) Mechanical state of the surrounding rock in a large-section clay tunnel.

A tunnel is excavated in a clay stratum, and the stress state and strain state of a point in the surrounding rock is the response of the internal forces and deformation of the surrounding rock under the generalized load. Therefore, the mechanical state of the surrounding rock can be characterized by its stress state and strain state.

The stress state and strain state of a point, M(x, y, z), in the surrounding rock of a clay tunnel can be expressed by the stress component and strain component of the element at that point (shown in Figure 2). σ = σ x τ x y τ x z τ y x σ y τ y z τ z x τ z y τ z ,   ε T = ( ε x ε y ε z γ x γ y γ z ) (2)

The following equation can express the relationship between the mechanical state of a point in the surrounding rock and its stress state and strain state. S = s   ( σ , ε ) (3)

S is the mechanical state tensor of a point in the surrounding rock after tunnel excavation. If the constitutive model is determined, the stress and strain state is the relevant variables. Then the mechanical state function of surrounding rocks can be expressed by one of the stress or strain states (for example, Eq. 4). S = s   ( σ ) o r S = s   ( ε ) (4)

For the clay tunnel surrounding rocks, failure of surrounding rock is caused by shear stress reaching the shear strength of clay. The stress state of a point in the clay surrounding rock can be replaced by the shear stress in Eq. 3, and the strain state is retained. The mechanical state function of the surrounding rock can be expressed by the following. S = s ( τ , ε ) (5)

Therefore, the actual mechanical state of a point in the surrounding rock under the unloading effect of tunnel excavation can be expressed as follows. S 0 = s ( τ 0 , ε 0 ) (6) Where τ 0 is the actual shear stress state at a point, and ε 0 is the actual strain state at the same point.

(3) The limit state of the surrounding rock in a large cross-section clay tunnel

The failure of the tunnel due to insufficient strength surrounding rock can be regarded as the mechanical state of the surrounding rock exceeding the limit state of bearing capacity. The failure of the tunnel due to excessive deformation of the surrounding rock can be regarded as the mechanical state of the surrounding rock exceeding the serviceability limit state. The tunnel limit state S _L has three situations of expression shown as follows.

(a) The stress state at a point in the surrounding rock reaches the limit, S L = s σ L ;

According to Eq. 1: when F s > 1 representing the surrounding rock in a stable state, the larger F s , the more stability reserves of the tunnel; when F s = 1 representing the surrounding rock in limit equilibrium, according to the stress-strain analysis can be seen at this time the surrounding rock in a state of yield; when F s < 1 represents the surrounding rock in a failure state.

Although Eq. 1 is a general expression, difficulties are encountered in the calculation process. The theoretical calculation of the mechanical state of the surrounding rock caused by tunnel excavation is relatively tricky. The shear strength can be simply determined by Coulomb’s law. Therefore, this paper uses the strength reduction method (Meng et al., 2019; Nie et al., 2019; Abra and Mahdi, 2022) and the strength-deformation parameter deterioration method (Huang et al., 2023) to calculate the tunnel stability reserve factor. The principle of the calculation of the tunnel stability reserve factor is simply represented by a Mohr circle and the strength envelope in Figure 3.

The factors of the stability of a large-section clay tunnel can be divided into internal and external factors. The tunnel stability reserve factor can be determined by Eq. 1. The mechanical conditions and the surrounding rock’s limit state also affect the tunnel stability. Equation 1 can be expressed by the following implicit function.

The complete set of influence factors of tunnel stability is U, and the groups of intrinsic and extrinsic factors are U i and U e , respectively ( U = U i ∪ U e ). Suppose t i ∈ U i and t e ∈ U e , then the binary implicit function of the tunnel stability reserve factor can be expressed as follows. F s = F t i , t e (7)

The internal and external variables in Eq. 7 contain at least one factor. The equation is not strictly a binary function.

(1) Internal factors affecting the stability of a large-section clay tunnel

The physical and mechanical properties of the surrounding rock (mainly including strength parameters, deformation parameters, weight, water content, et al.) are inherent and are internal factors affecting the stability of the tunnel. When the surrounding rock is clay, its physical and mechanical properties should conform to Coulomb’s law. And the shear strength of the surrounding rock can be determined by the strength parameters (Strength parameters are a collective term for the mechanical properties that affect the shear strength of clays) i.e., cohesion and angle of internal friction. The deformation of clay is determined by Young’s modulus and Poisson’s ratio via the generalized Hooke’s law. The deformation of the surrounding rock can be referred to as deformation parameters. The strength parameter and deformation parameter are the characterization of the inherent properties of materials.

Taking the limit state of the surrounding rock and the limit state of regular condition into account, the general expressions of the limit state of the surrounding rock at one point concerning the strength parameter and the deformation parameter are shown below. S L = s l c , φ , E , ν (8) Where c is the cohesion; φ is the angle of internal friction; E is Young’s modulus; and ν is Poisson’s ratio. c and φ are collectively referred to as the strength parameters; E and ν are collectively referred to as deformation parameters.

The above equations are supported by strength theory, and the constitutive relationship verifies their rationality. The clay’s bulk weight and water content may also be critical internal factors. To make the analysis more focused, their effects on strength and deformation are not considered here.

(2) External factors affecting the stability of a large cross-section clay tunnel

Tunnel excavation causes stress concentration or strain localization in the stratum in static equilibrium due to changes in structural form, boundary conditions, loading effects, et al. Moreover, the mechanical response (also called mechanical state) of the surrounding rock is the external condition and also is the key factor of the stability state of the surrounding rock. Therefore, the external factors affecting the stability of the tunnel are also the external conditions affecting the mechanical state of the surrounding rock, such as the action of additional loads on the surrounding rock, boundary conditions, the form of tunnel support, geometry, et al. These external conditions are closely related to the tunnel structure design plans, excavation plans, and support plans, which are collectively referred to as construction factors.

For a more focused subsequent analysis and to make the expressions more concise, the construction factors are here divided into three categories. For large cross-section clay tunnels, it is assumed that the set of factors related to the structural design scheme (e.g., burial depth, hole diameter, section shape, et al.) is A, which satisfies α ∈ A ; the set of factors related to the excavation scheme (e.g., excavation sequence, additional load, construction precipitation, et al.) is B, which satisfies β ∈ B ; the set of factors related to the support scheme (e.g., support form, support pressure, support distance, et al.) is C, which satisfies η ∈ C . The actual stress state of a point can be expressed as follows concerning Lame’s formula. τ 0 = τ 0 α , β , η (9)

The actual strain state at this point can be expressed as follows.

(a) When the stress-strain relationship is linear elastic

Regardless of the stress-strain relationship, Eq. 6 can be further expressed as: S 0 = s 0   c , φ , E , ν , α , β , η (12)

The joint analysis of Eq. 7 with Eqs. 8–12 leads to an implicit function in which more explicit variables can express the tunnel stability reserve factor. F s = f c , φ , E , ν , α , β , η (13)

When the surrounding rock material is the same, the stability state of the tunnel (which can also be considered as the actual mechanical state of the surrounding rock) must be different for different structural design plans, excavation plans, or support plans. Thus the tunnel stability index (i.e., the tunnel stability reserve factor) must be different. The engineers can estimate the current state characteristics of the tunnel via the quantitative evaluation of tunnel stability during the actual tunnel construction process. The stability evaluation guides the design, adjustment, or optimization of the construction can be planned to ensure the safety and stability of the tunnel, which is the most crucial goal pursued by engineers. Therefore, it is essential to determine the relationship between tunnel stability indicators and these main construction factors.

Large cross-sectional clay tunnels with little change in the distribution of the surrounding rock material can be regarded as interval sections of the same material. The tunnel stability reserve factor is closely related to the construction factors. At this point, Eq. 13 can be simplified to the following function. F s = f α , β , η (14)

Equation 14 reflects that we can use the tunnel stability reserve factor as a stability index while constructing the large cross-section clay tunnel. In developing a specific construction plan, this indicator can be used as a quantitative basis for the optimization and adjustment of the construction plan, and engineers can make construction decisions to optimize the safety and stability of the tunnel. The premise of all this is to determine the response of the tunnel stability reserve factor to these construction factors. α, β, and η represent one or more factors of a particular type. Therefore, to keep the discussion simple and general, α, β, η are still used here as a single variable for the analysis.

Any construction scenario for any large cross-sectional clay tunnels can be represented by a set of construction factors α i , β j , η k within the definition domain U. The tunnel stability reserve factor F _s for that set of construction scenarios can be defined as: F s = f α i , β j , η k (15)

According to the assumption of homogeneity of the surrounding rock material and physical-mechanical continuity, a higher-order partial derivative of F _s exists in the definition domain U. If the initial construction solution is denoted as α 0 , β 0 , η 0 , then the initial tunnel stability reserve factor can be denoted as F s 0 = f α 0 , β 0 , η 0 . The tunnel stability reserve factor F s can be denoted as the following equation according to Taylor’s theorem. F s = F s 0 + ∑ t = 1 n 1 t ! ( ∆ α ∂ ∂ α + ∆ β ∂ ∂ β + ∆ η ∂ ∂ η ) t f α 0 , β 0 , η 0 + R n (16)

Is the Lagrange residual term that can be used to estimate the error if there exists ξ ∈ 0,1 . R n = 1 n + 1 ! Δ t ∂ ∂ t + Δ σ t ∂ ∂ σ t + Δ σ s ∂ ∂ σ s 0 n + 1 F α r , β r , η r (17) Where, Δ α = α − α 0 Δ β = β − β 0 Δ η = η − η 0 , α r = σ + ξ Δ α β r = β + ξ Δ β η r = η + ξ Δ η (18)

Generally, as n increases, the higher the accuracy of the formula, but this will cause a rapid increase in the calculation and analysis of the computational effort, making its specific application difficult to promote; numerical analysis also shows that not higher-order Taylor formula will be able to achieve better accuracy, such as the Runger phenomenon. For the entire definition domain, complex functions are difficult to describe accurately with a simple Taylor formula. The first-order or second-order Taylor formula can follow the accuracy requirements in a relatively small sub-domain. The better fitting formula between the tunnel stability reserve factor and the many construction factors can be processed based on the above analysis. Equation 16 has a more apparent guiding meaning and practical value, which is further described as follows.

The initial construction scenario characterized by the main construction factors is α 0 , β 0 , η 0 , and the corresponding tunnel stability reserve factor F s 0 is characterized by the main construction factors, when the location of the large section tunnel is determined (i.e., the stratum is determined), the initial construction scenario. The stability of the tunnel under the initial construction scenario is thus quantitatively evaluated.

(a) If F s 0 ≤ 1 , then it indicates that the tunnel is difficult to maintain stability under the initial construction scenario.

To ensure a stable and safe tunnel construction, this faced two kinds of construction plan adjustments.

(a) When F s 0 ≤ 1 , the engineer needs to adjust the construction plan that α 0 , β 0 , η 0 → α i , β j , η k , so that the new construction plan at least meets F s > 1 .

Although we obtained the multivariate Taylor formula for tunnel stability index and construction factors, the procedure is still unknown. It is necessary to find their exact or approximate relationship equations.

Suppose the variables x 1 , x 2 , ⋯ , x n corresponding to any group of construction factors α , β , η are fitted, then the first-order response model of tunnel stability reserve factor F _s is shown in the following equation. F ^ s = ω 0 + ∑ i = 1 n ω i x i + δ (19)

The region corresponding to these construction factors is designed and calculated by the second-order model when the response surface of the tunnel stability reserve factor is close to or “curved”. In this case, there is a non-linear relationship between the stability reserve factor and these factors. In most cases, a second-order model is appropriate, as shown below. F ^ s = ω 0 + ∑ i — 1 n ω i x i + ∑ i — 1 n ω i i x i 2 + ∑ i < j ω i j x i x j + δ (22)

Additional experimental points are usually required in this region U 2 to achieve a better approximation. The estimation and optimization of the second-order response model also depend on a reliable experimental design scheme. Although there are other design methods, this paper uses the central composite design (CCD) as an example. Since the central composite design is robust to the assumption of strict linear effects, it also makes the RSM robust. The central composite design is a full or partial factorial design (usually with factors greater than five) with the addition of axis points n _a and center points n _c . The factor points are situated in the cube’s vertices, the centroids are located at the cube’s center, if the experimental design is represented as a ‘cube’, and the axis points are located on the cube’s axes in Figure 5.

The construction factors are coded and converted using Eq. 21. For CCD, each construction factor is now at three levels (i.e., high, medium, and low levels) of −1, 0, and 1. The distance of the axis point from the center point is assumed to be coded as λ . According to the center combination design, to satisfy the rotatability condition, the number of center points n _c needs to be appropriately selected to ensure that it is the predicted value with consistent uniform progress throughout the test area. Table 1 lists the number of experimental points of CCD with no more than five factors.

The number of experimental points in the CCD is more obviously influenced by some factors (as shown in Table 1). The distance λ of the axis point from the center point needs to be determined according to the rotatability, which can generally be calculated using the following equation. λ = 2 N f 4 (23)

N _f is the number of factors in the CCD.

A new tunnel somewhere on the Lanzhou Metro Line No. 3 multiplex is used as the study object. The tunnel is 970 m long and will cross the existing high-speed railway line at DK297 + 350 DK298 + 320. The tunnel axis is approximately 90° orthogonal to the high-speed railway line axis, and the minimum distance between the tunnel vault and the high-speed railway line is only 12 m. The new tunnel is mainly in the loess stratum. The tunnel in this scope is constructed by the shield method with a whole cross-section, and the excavation area is about 180 m², which is typical of a large cross-section shield tunnel.

The tunnel construction adopts the joint overrun support of a dense row of large pipe roofs and reinforcement bolts to ensure the stability of the tunnel construction and the minimum impact on the high-speed railway line. The row of large pipe roofs is constructed within a 150° angle and a 25 mm diameter reinforcement cage is inserted and filled with cement mortar to increase the stiffness (10 m in length, 180 mm in diameter, 10 mm in wall thickness, and 22 mm in diameter, 8 m in length and 1.2 m in spacing, in the plum-shaped arrangement).

The construction factors affecting tunnel stability in this paper mainly include the length of the pipe roof, the support pressure on the face, and the additional load on the ground. These construction schemes involve more complex conditions and are not convenient to use directly, a simplified approach is adopted in this paper to deal with them. For example, the length of the pipe roof is expressed as distance l from the face, the support pressure on the face is described as p, and the additional load on the ground is denoted as q, as shown in Figure 6.

The surrounding rock is a homogeneous clay layer and follows the ideal elastic-plastic constituted relationship. The numerical model is long enough in the axial direction, and the face is not less than five times the diameter of the hole from the left to the right boundaries to weaken the influence of the distal boundary on the analysis. The bottom of the tunnel arch is much larger than three times the tunnel’s diameter from the bottom edge. The other geometric dimensions of the numerical model (the diameter C = 20 m and the buried depth D = 14 m) and the initial clay physical and mechanical parameters are shown in Figure 6; Table 2.

The ideal tunnel model and the initial construction parameters are kept unchanged. The influence of the three construction factors (the length of the pipe roof l, the face support pressure, and the additional load on the ground as shown in Figure 6) on the stability of the face is analyzed. The region U ₁, U 1 = ( l , p , q l ∈ 0,20 , p ∈ 30,100 , q ∈ 0,60 is taken into account, where the unit of l is m, the unit of p is kPa, and the unit of q is kN/m.

The response surface method is used for numerical experimental design and function fitting. The finite element strength reduction method solves the tunnel stability reserve factor F _s . The finite element failure criteria are based on the unified use of the plastic zone penetration as the tunnel instability characteristics to avoid errors in the results. The different center point data acquisition is mainly achieved through variable grid size.

The full factorial experiments fitting the first-order response model were first arranged, and the obtained tunnel stability reserve factors are listed in Table 3. Note also that the tunnel stability reserve factor is calculated by the strength reduction method.

The resulting fit yields a first-order response model as: F s = 1.031 + 0.0505 l + 0.001 p − 0.006 q (24)

The fitted first-order response model, represented by the canonical variables, is F s = 1.37 + 0.505 ζ 1 + 0.015 ζ 2 − 0.18 ζ 3 (25)

Analysis of variance (ANOVA, R ² = 96.3%) shows that the first-order response model fits significantly in region U ₁. The ratio of the main effects of the pipe roof overrun support length, and support pressure on the face. The additional ground load will be further analyzed in the “curved” region U ₂ using the steepest ascent method as follows: b l : b p : b q = 1 : 0.029 : − 0.356 (26)

The steepest ascent method is conducted by a path of integer multiples of the step length Δ with the center point as the initial point, i.e., l ¯ = 10.0m, p ¯ = 15.0kPa, q ¯ = 30 kN/m.

Figure 7 shows the stability factor of safety for the face at each step along the steepest ascent path. The rate of increase in response is observed to be faster up to step 4. The rate of growth in the stability factor of safety remains essentially the same for each subsequent step. There is a significant “curved” after step 4, which requires a central composite design in the corresponding region U ₂. The region U ₂ for further fitting of the second-order response is U 2 = ( l , p , q l ∈ 18,26 , p ∈ 15.348,15.464 , q ∈ 24.30,27.15 (27)

Second-order response model for the tunnel stability reserve factor in region U ₂. F s = − 9.343 × 10 3 − 1.079 l + 1.209 × 10 3 p + 3.419 q + 1.08 × 10 − 1 l · p − 1.1 × 10 − 3 l · q + 4.785 × 10 − 2 p · q − 1.263 × 10 − 2 l 2 − 3.936 × 10 p 2 − 8.039 × 10 − 2 q 2 (28)

The results of model fit via the analysis of variance (ANOVA, R ² = 96.0%) show that the second-order model (Eq. 28) is a suitable approximation to the entire surface and a more accurate fit of the function F s = F l , p , q for the stability of the palm surface in the region U ₂.

Table 4 presents the results of the tunnel stability reserve factor of experimental points by CCD, and is used for second-order response fitting. The best-fitting expression can be obtained by effective tracing analysis of construction factors on these response surfaces (shown in Figure 8). The best-fitted model can be obtained in this numerical test, i.e., Eq. 28, which also reflects the strength of the interaction effect from the shape of contour lines in Figure 8. The interaction between the support pressure of the face and the additional load on the ground is more significant, which is consistent with the characteristics of their regression coefficients in Eq. 28. The contour line in Figure 8C is closer to the elliptical type, so the interaction between the support pressure of the face and the additional surface load is more significant. The response surfaces from Figure 8A ∼ (c) are a reasonable verification of the consistency and validity of Eq. 28 to the numerical tests. The corresponding construction plan can be determined more quantitatively according to this method in the early design and construction evaluation periods, thus ensuring the tunnel’s stability and safety.