Abstract
Due to the complexity of tunnels, accurate prediction of their loosening pressures in layered jointed rock strata is a very difficult engineering task. To recognize loosening patterns and estimate loosening pressures, numerical tests were employed in areas where tunnels were excavated in layered jointed rock strata. A total of 12 influential factors, including joints, tunnel depth, and strata, were considered in each of the numerical models. Three loosening patterns were found in the numerical testing: a ringent trumpet-shaped boundary, an arch-shaped boundary, and a closing-trend trumpet-shaped loosening zone. Empirical expressions for the loosening zone boundaries were further established and, in the form of the 12 influential factors, considered in the simulation. Given the boundary function, the loosening pressures were further deduced, which were categorized according to loosening pattern, i.e., ringent trumpet shape or arch shape, and the excavation condition of whether or not the embedded depth was deeper than the soft layer. Two case studies were used to test this method. The newly-proposed method was found to perform better than existing methods, with loosening pressure values that were slightly larger than, but very close to, actual measured field data.
Keywords
Highlights
Tunnels in layered jointed rock strata.
The loosening zone pattern was classified into three kinds.
Quantitative analysis on the weight of the influencing factors was analysed.
Predictions of the loosening zone boundaries were obtained.
Calculation formulae of loosening pressures were analytically solved.
Introduction
With the development of new infrastructure projects and buildings, as well as the renovation of existing facilities, the number and scale of tunnels constructed in jointed rock strata is rapidly increasing. Although rock itself may be homogenous, when we consider the rock mass over which we plan our construction, it may behave in an altogether different manner due to defects such as jointing, bedding planes, fissures, cavities, and other discontinuities. It is well known that the properties of jointed rock masses are very different from those of soils. Rock masses are discontinuous media with fissures, fractures, joints, bedding planes, and faults. These discontinuities may exist with or without gouge material and can cause significant or primary influence on rock mass strength. The spacing, orientation, roughness, and other characteristics of joints all have significant importance in terms of stability. Reliable characterization of the strength and deformation behaviour of jointed rocks is very important for the safe design of tunnels.
Nearly all of the existing methods for assessing earth/rock pressure, however, are based primarily on previous experience with sand or clay. Whether the support pressure calculations reflect the actual characteristics of the jointed strata is a key problem in tunnel design, and it is well understood that it may not be appropriate to apply earth pressure to the design of support systems in jointed rock masses.
Wang et al. 1 examined 8 factors affecting the loosening zone, including thickness of the weathered overburden and joint conditions. The failure modes of the loosening zone were classified and the formulae for the loosening pressures in shallow tunnels, which share the caving collapse loosening zone in layered jointed rock strata, were deduced using the upper bound theorem. The proposed calculation formulae were restricted, however, by the specific structure size and strata properties of the Xinggongjie metro station in Dalian City, China. In addition, the long calculation procedure of that method is more complicated than that of traditional methods.
Focusing on a specific shallow tunnel in layered jointed rock strata, and employing the orthogonal array testing strategy (OATS) and the distinct element method (DEM), this study extends Wang’s work 1 by considering 12 influential factors, such as the structural sizes of tunnels, configurations of the rock joints, thickness of the overburden, properties of the joints and rock strata, and tunnel depth. First, a total of 384 numerical models were established using the universal distinct element code (UDEC). 2 Through comprehensive analysis of the test results, the failure modes of the loosening zone were classified. Second, the data envelopment analysis (DEA) 3 theory was introduced and used to analyse the weights or influence degree of the influential factors. Based on this analysis, the major influential factors were chosen and used in the one-dimensional nonlinear and multiple linear regression models. Subsequently, the prediction formulae for the boundaries of the loosening zone were obtained. Third, analytical models were developed using the stress transformation theory, and the formulae for the loosening pressures were established. Finally, 2 case studies were performed in order to verify the proposed method.
Methods
Given the expansive scale of contemporary underground engineering projects, the stability problems of tunnels have become extremely important. In comparison to old tunnels, proposed tunnels in China and abroad are longer, wider, and have heavier overburdens. Therefore, it is necessary to clearly understand the action mechanisms of the factors affecting the loosening zone of the surrounding rock. The influential factors can be subdivided into 3 types: (1) geological factors such as the physical and mechanical properties of the rock strata; (2) factors such as the scale, depth, and height-span ratio; (3) execution factors such as construction methods and support measures. It is well known that geological formations and rock structures govern the surrounding rock stability. Until now, many studies have examined the influential factors of tunnel stability, such as those found in Liu and Song, 4 Xu and Wang, 5 and Wang and Zhu.6,7 These studies have the following limitations: (1) In general, most of the research focused on a specific case, and only a single factor—such as tunnel depth—was analysed, as opposed to conducting a multifactorial analysis; (2) the analysis of influence degree was qualitative, when it is necessary to explore quantitative measurements of the influential factors; (3) these studies are mostly applicable to soil, while research needs to be conducted for jointed rock strata; (4) factors such as joint state and overburden thickness were not taken into account.
In this study, the excavation of a tunnel is simulated using the DEM. In order to quantitatively evaluate the effect of upper-level factors on tunnel construction, a total 384 numerical tests were carried out. Based on the results from the numerical simulation, the major loosening zone patterns could then be recognized and the weight of each factor solved using DEA theory. Further, the shape, range, and expressions of the loosening zone boundaries could be empirically linked to these factors.
Numerical simulation
Numerical simulation method
Numerical simulation has become a useful tool in the design of underground excavations, in which the magnitude and distribution of earth pressure on support systems can be investigated under various constraints such as limited amounts of time, cost, and space.
Useing remeshing, the finite element method (FEM) is the most widely applied technique for stratum settlement simulations. 8 The assumption of the FEM is that the rock mass is a continuous medium which ideally can be represented by an equivalent continuum constitutive model. However, since rock masses inherently contain discontinuities of many types such as fissures, fractures, joints, bedding planes, and faults, improvements are needed to simulate the multiple growths of discontinuities. Boundary element methods (BEMs) have been applied in both continuum formulations and discontinuous formulations. However, BEMs are difficult to use when solving nonlinear problems, which involve integrating the nonlinear terms with singularities. 1 Other attempts have been made using meshless methods, although this approach has yet to attain broad application.9,10
Since the discontinuous method can simulate the discontinuous nature of the rock strata failure process, we believe this method is the most logical choice. Cundall 11 and Cundall and Hart 12 developed the DEM to model the responses of discontinuous rock masses. This technique can directly approximate the blocky structure of jointed rock using arbitrary polyhedrons and can recognize and handle new contacts automatically during the calculation. At present, the DEM is rapidly gaining in popularity for the behavioural analysis of tunnels excavated in jointed rock masses.13–15
In this study, the UDEC based on the DEM was adopted to examine the excavation process of tunnels in jointed rock. The UDEC allows large displacements between blocks, can display the displacement and velocity of rigid or deformable blocks, and more. It also provides high-resolution graphics of the deformation process.
During the simulation, block systems are considered in which the Mohr-Coulomb failure criterion is employed for the rock block and the Coulomb slip model is implemented for the rock joints, i.e., sliding contacts among blocks. In most cases, the calibration of the parameters involved in the Coulomb model is much easier than the calibration of other complex models, e.g., Barton-Bandis joint models.
16
Complex joint models are empirical expressions requiring more in-depth knowledge of joint behaviour, which may not be suitable for general numerical parametric analyses. The Coulomb slip model can be expressed as follows:
Referring to Goodman,
17
considering a jointed rock mass with uniform joint spacing smean, kn and ks are given by:
When the joint shear stress exceeds the joint shear strength, which is a combination of cohesive (cJ) and frictional (φJ) strength, the joint loses strength, causing the displacement of the block. The joint model is the most suitable approach for parametric analyses in general engineering and can approximate the behaviour of joints.
Numerical model description and testing programme
Most published studies of stability problems of jointed rock masses only consider 2 joint sets, e.g., Yeung and Leong, 18 Solak and Schubert, 19 and Delisio et al. 20 Based on the Chinese rock mass subclassification system,21–25 4 independent sets of joints were considered, including 2 sets of inclined joints with different orientations for both soft and hard rock, 1 set of horizontal joints for soft rock, and another set of horizontal joints for hard rock.
Most published studies exploring the stability problems of jointed rock masses have only considered 2 joint sets.18–20 Based on the Chinese rock mass subclassification system,21–25 the qualitative description of the subclassification of surrounding rock mass grades IV and V (Table 1) can be obtained, from which the number of joint sets and the range of joint spacing can be determined. In this study, the rock masses were divided into 2 layers, i.e., the upper soft stratum and the lower hard stratum, as shown in Figure 1. For the soft and hard rock strata illustrated in Figure 1, class IV rock including the IV1, IV2, IV3 subclassifications, and class V rock, including the V1, V2 sub-classifications, were considered. Six stratum combinations were given as follows: V1–IV1, V2–IV1, V1–IV2, V2–IV2, V1–IV3, and V2–IV3. The properties of the rock strata are listed in Tables 2 and 3.

Schematic diagram of joints.
Qualitative subclassification methods for surrounding rock mass grades IV and V.
Physical and mechanical parameters of the rock masses and joints.
Note: Ee and ve are the Young's modulus and Poisson's ratio for the rock mass, respectively; φJ and cJ are the friction angle and cohesion of the joints, respectively.
Physical and mechanical parameters of the rock.
Note: γr, Er, vr, φr, and cr, are the density, Young's modulus, Poisson's ratio, friction angle, and cohesion of the intact rock material, respectively.
The ranges of the other factors were determined according to the subclassification system, engineering experiences, field investigations, and observations.26–32 The span of the tunnel, W, varied from 5–26 m; the height-span ratio of the tunnel, λ, varied from 0.5–2.0; the thickness of the soft overburden, h, ranged from 5–50 m; and the depth of the tunnel, H, varied from 10–100 m. α1 is the angle of the transfixion inclined joints (10° ≤ α1 ≤ 80°); α2 is the angle of the intermittent inclined joints ( − 80° ≤ α2 ≤ −10°); the lateral pressure coefficient, k, ranged from 0.5–2; s1 represents the joints spacing of the transfixion inclined joints and the intermittent inclined joints; s2 denotes the spacing of horizontal joints in the regolith stratum; and s3 denotes the spacing of horizontal joints in the bedrock stratum. The values of s1, s2, and s3 vary with stratum class, as shown in Table 4.
Range of joint spacing values.
The numerical tests were designed using OATS theory. 1 There were 12 factors: H, W, λ, h, k, α1, α2, s1, s2, s3, and the stratum parameters of the upper and lower strata, Ds and Dx, respectively. The range of each factor's value was divided into 4 equal parts, i.e., each factor had 4 levels. It is worth noting that each level of Ds and Dx was treated as a group during the tests, as shown in Tables 5–7, where Ds1, Ds2, Ds3, and Ds4 are the 4 levels of the properties of the upper soft stratum; Dx1, Dx2, Dx3, and Dx4 are the 4 levels of the properties of the lower hard stratum; r1, φ1, c1, E1, and v1, denote the density, friction angle, cohesion, elastic modulus, and Poisson's ratio of the upper soft stratum, respectively; r2, φ2, c2, E2, and v2 denote the density, friction angle, cohesion, elastic modulus, and Poisson's ratio of the lower hard stratum, respectively; cJ1 and φJ1 represent the joint cohesion and joint friction angle in the upper soft stratum, respectively; and cJ2 and φJ2 represent the joint cohesion and joint friction angle in the hard stratum, respectively. Taking the stratum combination V1–IV1 as an example, the range and level of the factors are listed in Table 8. For each stratum combination, 64 test scenarios were determined using the orthogonal array table L64(1212); 1 thus, there were 384 test scenarios in total for the 6 stratum combinations.
Physical and mechanical parameters of the rock and joints in the upper soft stratum.
Physical and mechanical parameters of the rock in the lower hard stratum.
Physical and mechanical parameters of the joints in the lower hard stratum.
Range and level of the factors for the stratum combination V1–IV1.
Classification of the loosening zone
The collapse area can be obtained after the numerical calculation. Boon et al. 33 examined the collapse patterns of tunnels in complex jointed rock masses consisting of 3 independent sets of joints, and reached the following conclusions: (1) It is possible for all of the joint sets to participate in the collapse mechanism and (2) the collapse mechanisms are numerous. Therefore, it emerges that the prediction of the exact failure location is unlikely to be successful in practice since the collapse location is highly sensitive to the joint patterns and there could be more than one failure mechanism occurring simultaneously. For these reasons, it appears as if the qualitative characterization of collapse failure patterns using specific descriptions is unattainable.
Since the strength and stiffness of rock blocks are large, damage to the surrounding rock is mainly due to the opening and slipping of the joints. Hence, the loosening area can be regarded as the slipping or opening zone. The UDEC can display the slipping and opening zones clearly, as shown in Figure 2. As can be seen in the figure, the slipping zone area is normally larger than the opening zone. Hence, the loosening zone can conservatively be regarded as the slipping zone.

Loosening zone patterns for tunnels of different depths.
Referring to Figure 2, the loosening zone of the surrounding rock masses can be classified into 3 patterns: the ringent trumpet-shaped loosening zone, the closing-trend trumpet-shaped loosening zone, and the arch-shaped loosening zone. The number of each kind of loosening zone via the 6 stratum combinations is listed in Table 9. This reveals that (1) the number of ringent trumpet-shaped loosening zones decreases with increasing tunnel depth, while contrary to this, the number of arch-shaped loosening zones increases with increasing tunnel depth, and the closing-trend trumpet-shaped loosening zone is the transition state between the other 2 kinds of loosening zone; (2) the number of ringent trumpet-shaped loosening zones is 217 (57 percent of total), and the number of arch-shaped loosening zones is 143 (37 percent of total); (3) the models mostly exhibit ringent trumpet-shaped loosening zones when the depth of the tunnel is less than or equal to 70 m, and models mostly exhibit arch-shaped loosening zones when the tunnel depth is less than or equal to 100 m or greater than 70 m. Therefore, this study focused on the ringent trumpet-shaped and arch-shaped loosening zones.
Number of the 3 types of loosening zones at different depths.
Factors influencing the loosening zone
Data envelopment analysis
DEA is a linear programming technique. 3 It is a non-parametric method used in the estimation of production functions and has been utilized extensively to estimate measures of technical efficiency in a range of industries. 34 Like the stochastic production frontiers, DEA estimates the maximum potential output for a given set of inputs, and has primarily been used in the estimation of efficiency. A key advantage of DEA over previous approaches is that it more easily accommodates both multiple inputs and multiple outputs. Some of the advantages of DEA are: (1) There is no need to explicitly specify a mathematical form for the production function; (2) it has proven to be useful in uncovering relationships that remain hidden from other methodologies; (3) it is capable of handling multiple inputs and outputs; (4) it is capable of being used with any input-output measurement; and (5) the sources of inefficiency can be analysed and quantified for every evaluated unit.
A range of DEA models have been developed that measure efficiency and capacity in different ways. These largely fall into the categories of being either input-oriented or output-oriented models. The CCR model is named after its developers Charnes, Cooper, and Rhodes. This was the first and fundamental DEA model, built on the notion of efficiency as defined by the classical engineering ratio. This model calculates an overall efficiency for the unit in which both its pure technical efficiency and scale efficiency are aggregated into a single value. The obtained efficiency is never absolute as it is always measured relative to the field.
As mentioned in the previous section, based on OAT theory, 384 numerical models were calculated, with 12 influential factors included in each. These influential factors, i.e., the input factors, were controllable; however, the results of the numerical tests—the outputs—such as the area and height of the loosening zone, were uncontrollable. Therefore, in the present study, based on the input-oriented CCR model, the weights of the individual influential factors were analysed.
Quantitative analysis of the weight of each influential factor
For both the ringent trumpet-shaped loosening zone and the arch-shaped boundary, the input factors included H, W, λ, h, k, α1, α2, s1, s2, s3, Ds, and Dx. The precise values of each of these factors are listed in Tables 5–8. As shown in Figure 3, X is defined as the lateral length of the loosening zone on the ground surface; L1 and L2 are the widths of the horizontal earth strip above tunnel vault; A1 and A2 are the areas of the loosening zone; and V1 is the height of the loosening zone above the vault of the tunnel. For the ringent trumpet-shaped loosening zone, the outputs were X, L1, and A1; for the arch-shaped loosening zone, the outputs were A2, L2, V1, and V2. The statistical data for X, L1, A1, A2, L2, V1, and V2 are listed in the supplementary material (Tables S1–S3).

Schematic diagram of model outputs.
The efficiency of each of the influential factors was analysed using MaxDea code; results of this analysis are listed in Tables 10 and 11. These results show that: (1) For the ringent trumpet-shaped loosening zone, the relative degree of influence of the factors was W > H > λ > α1 > h > Ds > s1 > k > s2 > Dx > α2 > s3; for the arch-shaped loosening zone, the relative degree of influence was W > λ > H > α1 > Dx > s1 > k > s3 > Ds > s2 > α2 > h. This is due to the fact that, with the exception of Dx and Ds, the other 10 factors all share the same value in the 6 stratum combinations; (2) for both types of loosening zone, the span of the tunnel, W, is the most important influential factor. The state of the joints also has significant influence that cannot be neglected in the failure of the loosening zone; (3) for the ringent trumpet-shaped loosening zone, factors such as α2 and s3, which are associated with the hard stratum, exert less influence on the failure; correspondingly, for the arch-shaped loosening zone, factors such as α2 and h, associated with the soft stratum, exert less influence.
Factor weights for the ringent trumpet-shaped loosening zone.
Factor weights for the arch-shaped loosening zone.
Prediction of loosening zone boundaries
It is well known that the range of the loosening zone has a significant impact on the construction safety and support of tunnels; hence it is of great importance to derive a formula for the loosening zone boundary, from which the loosening zone range can be defined. The previous section analysed the influential factors of the loosening zone, revealing that different stratum combinations have different factor weights. In this study, the factors with weights greater than 2% were considered to be important parameters that could not be ignored. On this basis, the one-dimensional nonlinear and multiple linear regression models were constructed, and the prediction formulae of the boundaries of the ringent trumpet-shaped and arch-shaped loosening zones were deduced. The first step was to identify the functions for the 2 types of loosening zone boundaries.
The functions of the loosening zone boundaries
Curve fitting is one of the most powerful and most widely used analytical tools in Origin code. Origin code provides tools for linear and nonlinear curve fitting, along with validation and goodness-of-fit tests. Let the expression y = f(x) indicate any element belonging to the set of curves which define a possible range of the loosening zone in the Cartesian reference frame x-y.
For the ringent trumpet-shaped loosening zone exhibited by 217 models, the right bottom corner of the tunnel is used as the origin of coordinate system, as shown in Figure 4.

Schematic diagram of models with ringent trumpet-shaped loosening zones.
The right boundary of the ringent trumpet-shaped loosening zone shown in Figure 4 can be simulated by 2 functions, the power function (y = axb) and the parabolic function (y = a + bx + cx2). For each test, the undetermined coefficients, i.e., a and b in the power function and a, b, and c in the parabolic function, can be determined with regression tools by utilizing Origin code. As examples, the boundaries of loosening zones from numerical tests No. 8, No. 20, No. 28, and No. 35 are shown in Figure 5; regression results are listed in Table 12. Due to its higher coefficient and lower number of parameters, the power function was chosen as the boundary function of the ringent trumpet-shaped loosening zone in this study.

Fitting comparison of the power function and the parabolic equation.
Fitting results for the power function and parabolic equation.
Note: a, b, and c are the values of the coefficients; σa, σb, and σc are the standard errors of the coefficients; R is the correlation coefficient of the equations.
Similarly, for the 143 models with the arch-shaped loosening zone, as shown in Figure 6, the expression of the boundary can be simulated by 3 functions: the inversely proportional function (y = 1/[a + bx]), the power function (y = a + bxc), and the parabolic function (y = a + bx + cx2). The boundaries of the loosening zones from numerical tests No. 22, No. 21, No. 53, and No. 59 are shown in Figure 7; regression results are listed in Table 13. Due to its higher coefficient and lower number of parameters, the inversely proportional function was chosen as the boundary function of the ringent trumpet-shaped loosening zone in this study.

Schematic diagram of models with arch-shaped loosening zones.

Fitting comparison of the inversely proportional function, the power function, and the parabolic equation.
Fitting results for the inversely proportional function, power function, and parabolic equation.
Note: a, b, and c are the values of the coefficients; σa, σb, and σc are the standard errors of the coefficients; R is the correlation coefficient of the equations.
Prediction formulae of the loosening zone boundaries
Prediction formulae of typical geometric quantities
For the 217 models with the ringent trumpet-shaped loosening zone, according to Tables 5–8, one-dimensional nonlinear regression models of X (the lateral length of the loosening zone on the ground surface) and L1 (the width of the horizontal strip of earth above the tunnel vault) can be established using the least squares method. As an example, for the V1–IV1 stratum combination, the one-dimensional nonlinear regression formulae of X are listed in Table 14, where 10 m ≤ H ≤ 40 m; 5 m ≤ W ≤ 26 m; 0.5 ≤ λ ≤ 2; 5 m ≤ h ≤ 50 m; 0.5 ≤ k ≤ 2; and 10° ≤ α1 ≤ 80°. For the stratum combinations V1–IV1, V1–IV2, and V1–IV3, the range of s1 is 0.4 m ≤ s1 ≤ 1.0 m; for the stratum combinations V2–IV1, V2–IV2, and V2–IV3, the range of s2 is 0.2 m ≤ s2 ≤ 0.4 m. When the values of the stratum parameters are approximate to Ds1 (see Table 5), then Ds takes the value 1; analogously, Ds and Dx can take the values 1, 2, 3, or 4, corresponding to different stratum parameters. Statistically speaking, the values of the correlation coefficient, R, are all greater than 0.900, indicating that all of the empirical formulae are correct.
One-dimensional nonlinear regression formulae of X for the stratum combination V1–IV
Suppose X and 1.005
H
, W, λ2, λ, h, e0.222k, tan2α1, tanα1, s12, s1, s22, s2, Ds, and Dx have a linear relationship. Through the construction of multiple linear regression models, the value of each undetermined coefficient can be calculated using the least squares method and the calculation formula of X obtained:
For other stratum combinations, the multiple linear regression equations of X and L1 are listed in Tables 15 and 16. Similarly, the multiple linear regression equations of L2, V1, and V2 can be obtained (Tables 17–19). The F distribution values of these equations are listed in Tables 15–19; the rationality of the above equations was also verified using ANOVA.
Multiple linear regression equations of X.
F*: F distribution values.
Multiple linear regression equations of L
F*: F distribution values.
Multiple linear regression equations of L
F*: F distribution values.
Multiple linear regression equations of V
F*: F distribution values.
Multiple linear regression equations of V
F*: F distribution values.
Prediction formulae
In the case of the ringent trumpet-shaped loosening zone shown in Figure 4, the boundary function of the loosening zone was considered to be y = axb. Once the prediction formulae of X and L1 had been derived, the coefficients, a and b, could be calculated:
Calculation formulae for loosening pressures
Once the boundary of the loosening zone was obtained, the force field could be determined using the principle of stress transformation. Let r1 represent the density of the soft stratum, cJ1 denote the joint cohesion of the soft stratum, φJ1 represent the joint internal friction angle of the soft stratum, r2 represent the density of the hard stratum, cJ2 denote the joint cohesion of the hard stratum, φJ2 represent the joint internal friction angle of the hard stratum, h denote the depth of the soft stratum, H denote the depth of the tunnel, and W represent the span of the tunnel. The analytical formulae of the loosening pressures were deduced for the following 2 conditions: (1) embedded depth deeper than the soft layer (H > h); (2) embedded depth shallower than the soft layer (H ≤ h).
Condition 1: The embedded depth is deeper than the soft layer (H > h)
Loosening pressures of models with ringent trumpet-shaped loosening zones
As shown in Figures 8 and 9, the straight line, M1N1, denotes the interface between the soft and hard strata, which divides the loosening zone above the tunnel vault into 2 parts: the geometric regions E1M1N1B1 and M1R1K1N1.

Calculation model (No. 1).

Calculation model (No. 2).
The geometric region D1U1V1C1 denotes the differential stratum strip. Force field analyses were carried out; details of the derivation of the support pressures are illustrated in Appendix 1. The computational formulae of the vertical pressure of the tunnel, q, could then be obtained:
For the calculation of the lateral pressure, according to traditional approaches such as Terzaghi theory, 35 the calculation formulae are given on the basis of Rankine's ratio.
However, it is well known that Rankine's ratio (the active earth pressure coefficient) is applied to earth pressure with the conditions: (1) The smooth surface which is the interface between the wall and soil is frictionless; (2) the supported soil is homogeneous and isotropic. In addition, when the critical slip planes are oriented at 45° ± φ/2, the active Rankine's ratio is ka = tan2(45°-φ/2). Krynine, 36 Ladanyi and Hoyaux, 37 and Handy 38 found that since ka was derived from the assumption that the horizontal and vertical stresses are the principal stresses, using ka is only valid when there is no friction (shear stress) between the till material and the sides of the ditch. However, Terzaghi theory posits that friction exists at the sides of the ditch. Therefore, the active Rankine's ratio, ka, cannot be used. Wang and Zhu 39 also verified this perspective.
In order to avoid the unrealistic assumption of Rankine's active earth pressure, the slip angle is corrected using the following formula (referring to Figure 8):
Loosening pressures of models with arch-shaped loosening zones
As is shown in Figures 10 and 11, the straight line, M2N2, denotes the interface of the soft and hard strata, while the geometric region, D3U3V3C3, denotes the differential stratum strip. Force field analyses were carried out; details of the derivation of the support pressures are illustrated in Appendix 2.

Calculation model (No. 3).

Calculation model (No. 4).
The vertical force, q*, can be written as
The lateral pressure can be written as
Condition 2: The embedded depth is shallower than the soft layer (H < h)
Loosening pressures of models with ringent trumpet-shaped loosening zones
As shown in Figure 12, when h ≥ H, the tunnels are in homogeneous rock masses. The derivation procedure for q is the same as that used for q0, which was described in section 4.1.1. The calculation formulae of the loosening pressure can be directly written as:

Calculation model (No. 5).
Loosening pressures of models with arch-shaped loosening zones
As shown in Figure 13, when h ≥ H, the tunnels are in homogeneous rock masses. The derivation procedure for q* is the same as that used for q0*, which was described in section 4.1.2. The calculation formulae of the loosening pressures can be directly written as:

Calculation model (No. 6).
Examples and discussion
Case studies
In many regions of China, rock strata are frequently present at excavation sites, since the soil layer is relatively shallow and the construction of high-rise buildings or underground facilities requires larger and deeper excavations. In this section, 2 actual tunnel projects were utilized as case studies. In addition to the proposed method, for models with ringent trumpet-shaped loosening zones, Bierbäumer's theory 40 and Terzaghi's theory 1 were used; for models with arch-shaped loosening zone, Balla's theory 41 and Fraldi and Guarracino's theory 1 were utilized.
Case 1
The No. 2 Dalian metro line of Dalian, China is 42.56 km long and has 29 underground stations. Xinggong Street is one of these underground metro stations. This station is embedded in heavily jointed rock strata. The lithology of the area consists primarily of a succession of strongly weathered calcareous slate (soft stratum) and moderately weathered calcareous slate (hard stratum). The mechanical parameters of the rock mass, joints, and rock are listed in Table 20.
Physical and mechanical parameters of the rock mass and joints.
SCS*: strongly weathered calcareous slate; MCS*: moderately weathered calcareous slate; γe, Ee, ve, φe, and ce are the density, Young's modulus, Poisson's ratio, friction angle, and cohesion of the rock mass, respectively; φJ and cJ are the friction angle and cohesion of the joints, respectively; γr, Er, vr, φr, cr, and pr are the density, Young's modulus, Poisson's ratio, friction angle, cohesion, and uniaxial compressive strength of the intact rock material, respectively.
As shown in Figure 14, this tunnel is 16.865 m in height (T) and 21.6 m in width (W), with a cross section greater than 280 m2. From a field survey, the joint spacing (s1) and angle (α1) of the transfixion inclined joints were determined to be 0.8 m and 55°, respectively. The joint spacing of the horizontal joints (s2) in the upper soft stratum was 0.35 m. The lateral pressure coefficients of the soft stratum (k1) and the hard stratum (k2) were 0.65 and 0.70, respectively. The thickness of the strongly weathered calcareous slate (h) was 21 m.

Tunnel geometry of the Xinggong street metro station in Dilian, China.
It is worth pointing out that since the strata layers extending from the ground surface to the bottom of the tunnel are soft in the upper regions and hard in the lower regions, the calculation should take into account the characteristics of the layered strata. Thus, the values of the rock density, cohesion, and internal friction were weighted according to the thickness of each layer when using Bierbäumer's theory and Terzaghi's theory.
Letting H take different values ranging from 24–40 m, the vertical pressure (q) calculated using the two theories can be obtained (Figure 15).

Results comparison No. 1 (calculated using Bierbäumer's theory, Terzaghi's theory, and the proposed method for shallow tunnels in heterogeneous rock masses).
As mentioned in the previous section (see Table 9), when H ≤ 40 m, 96 percent of the models had a ringent trumpet-shaped loosening zone; therefore, the loosening zone in this case study was considered to be ringent trumpet-shaped. Utilizing Eqs. (11) and (12), the values of q calculated using the proposed method could also be obtained (see Figure 15)
From Figure 15, it can be concluded that the pressure values calculated using the proposed method were much smaller than those calculated using Bierbäumer's theory and Terzaghi's theory. In addition, when H = 25 m, the field monitoring data indicated that q = 180 kPa, which was in reasonable agreement with the values calculated using the proposed method (in Figure 15, when H = 24 m, q = 177.34 kPa; when H = 26 m, q = 182.54 kPa).
For tunnels in homogeneous rock masses, the Xinggong Street metro station could still be used as an example, with the condition that H = h, while all other parameters remained unchanged. Similar conclusions could also be obtained (see Figure 16) by virtue of Eqs. (20) and (21). Therefore, it can be concluded that, compared to the empirical method, the method proposed in this study is a valid and efficient approach for the determination of loosening pressures.

Results comparison No. 2 (calculated using Bierbäumer's theory, Terzaghi's theory, and the proposed method for shallow tunnels in homogeneous rock masses).
Case 2
Zhenlong tunnel is a railway tunnel located in Zunyi City, China. The overall length of this tunnel is 1573 m, and its depth varies from 75–95 m. The lithology of the strata primarily consists of a succession of strongly weathered limestone (soft stratum) and moderately weathered mudstone (hard stratum), in which the thickness of the soft stratum is 20 m. The mechanical parameters of the rock mass, joints, and rock are listed in Table 21. The meanings of the symbols are identical to those described in Table 19.
Physical and mechanical parameters of the rock mass and joints.
As shown in Figure 17, this tunnel is 8.13 m in height (T) and 5.7 m in width (W). From a field survey, the joint spacing (s1) and angle (α1) of the transfixion inclined joints were determined to be 0.5 m and 61°, respectively. The joint spacing of the horizontal joints (s2) in the upper soft stratum was 0.31 m. The lateral pressure coefficients of the soft stratum (k1) and the hard stratum (k2) were 0.70 and 1.0, respectively. The thickness of the strongly weathered calcareous slate (h) was 18 m.

Zhenlong tunnel geometry.
Letting H take different values ranging from 75–95 m, the vertical pressure (q) can be calculated using Balla's theory. It is worth pointing out that since the strata layers extending from the ground surface to the bottom of the tunnel are soft in the upper regions and hard in the lower regions, the calculation should take into account the characteristics of the layered strata. Thus, the values of fH, fw, fc, and γ were weighted according to the thickness of each layer when using Balla's theory. For Fraldi and Guarracino's method, the calculation should also take into account the characteristics of the layered strata. As mentioned in the previous section (see Table 9), when H ≥ 70 m, 83 percent of the models had an arch-shaped loosening zone; therefore, the loosening zone in this case study was considered to be arch-shaped. Utilizing Eqs. (16) and (17), the values of q calculated using the proposed method could also be obtained (see Figure 18).

Results comparison No. 2 (calculated using Balla's theory, Fraldi and Guarracino's method, and the proposed method for tunnels in heterogeneous rock masses).
A comparative analysis indicated that: (1) The pressure values calculated using the proposed method were smaller than those calculated using either Balla's theory or Fraldi and Guarracino's method. In addition, the differences among the values calculated using the 3 methods increased with tunnel depth; (2) the values calculated using the proposed method were greater than the field monitoring data, and all of the relative differences were less than 5% (Table 22 and Figure 14).
Relative differences between the results calculated using the proposed method and the field monitoring data.
R*: Relative difference.
For tunnels in homogeneous rock masses, the Zhenlong tunnel could still be used as an example, with the condition that H = h, while all other parameters remained unchanged. Similar conclusions could then also be obtained (see Figure 19). Therefore, it can be concluded that the method proposed in this study is a valid and efficient approach for the determination of loosening pressures.

Results comparison No. 4 (calculated using Balla's theory, Fraldi and Guarracino's method, and the proposed method for tunnels in homogeneous rock masses).
Discussion
It is necessary to discuss the mechanism and result differences between the proposed method and the traditional methods.
First, the existing traditional methods for assessing earth/rock pressure are primarily based on previous experience with sand or clay soils. Coefficients used in these methods may include the conditions (such as tunnel depth) and the properties (such as friction angle) of the geomaterial, but for layered jointed rock strata, the loosening pressures of a tunnel are obviously affected by other factors such as the configuration of the rock joints. Given the different characteristics of rock joints and the variability of strata, applying these methods would introduce errors into the results.
Second, there are some inherent uncertainties in these traditional methods. In particular, some of the assumptions they employ may not be realistic under certain circumstances. For example, in Terzaghi's theory, the sliding surfaces are assumed to be vertical and the cross section of the loosening zone of the surrounding rock is depicted as a rectangle. Thus, the shape and boundary of the loosening zone are only affected by tunnel depth, tunnel height, and internal friction angle of the stratum. In Bella's theory, many assumptions based on the researcher's experience were made. The defect in this theory is that both the loosening patterns and pressures of the lateral rock mass adjacent to the tunnel are neglected since the sliding surfaces begin from a corner point. Therefore, it is critical for researchers to modify the limitations of the assumptions.
Third, the benefits of using traditional methods include speed and ease of use, provided that engineers select a method that was derived from ground conditions and construction techniques similar to those of their proposed project. If not, errors will occur. 41 Rock masses are discontinuous media, and the properties of the jointed rock masses are very different from those of soils. Given that most of the traditional methods have only been applied to soil, these methods should be applied in conjunction with other analytical approaches, as well as field monitoring, in order to confirm their results.
Conclusions
Loosening patterns and pressures are of great concern in tunnel construction. In order to recognize the loosening pattern and estimate the loosening pressure, a comprehensive numerical analysis was performed utilizing 2D UDEC code. In the numerical testing, tunnel excavations in layered rock strata were simulated, in which a lower hard-jointed rock layer was overlaid with an upper soft rock layer. A total of 12 factors, including the depth of the tunnel (H), the thickness of the soft stratum (h), the span and the height-span ratio of the tunnel (W and λ, respectively), the lateral pressure coefficient (k), the spacing of 2 sets horizontal joints in the upper soft stratum and lower hard stratum (s2 and s3, respectively), the angle of the transfixion inclined joints (α1), the angle of the intermittent inclined joints (α2), the spacing of the 2 sets inclined joints (s1), and the stratum parameters of the upper and lower strata (Ds and Dx, respectively) were considered in the numerical tests. A total of 384 numerical models were conducted. The results indicate that:
Three loosening patterns were obtained: a ringent trumpet-shaped loosening zone, a closing-trend trumpet-shaped loosening zone, and an arch-shaped loosening zone. The predominance of a particular loosening pattern was discovered to be a function of both the embedded tunnel depth and the dip angle of the joints. Most tunnels (96%) with embedded depths less than 40 m exhibited ringent trumpet-shaped loosening zones, while most tunnels (83%) with embedded depths ranging from 70–100 m exhibited arch-shaped loosening zones. Among the 12 factors considered, W, H, and λ were the most important. Specifically, for ringent trumpet-shaped loosening zones the relative degrees of influence of the factors were W > H > λ > α1 > h > Ds > s1 > k > s2 > D x > α2 > s3, while for arch-shaped loosening zones the relative degrees of influence were W > λ > H > α1 > Dx > s1 > k > s3 > Ds > s2 > α2 > h. Regression analysis conducted on the boundary shapes of the loosening zones resulted in 2 empirical expressions: one for ringent trumpet-shaped loosening zones and another for arch-shaped loosening zones. These empirical expressions can be used to predict the shapes of loosening zones for tunnels constructed in layered rock strata. Given the shape of a loosening zone, its loosening pressure can be analytically solved. Four series of formulae were deduced according to the loosening pattern, i.e., ringent trumpet-shaped or arch shaped, and an excavation condition in which the embedded depth was either deeper or shallower than the soft layer.
The newly proposed method was applied to two practical tunnel projects. The loosening pressures predicted by the newly proposed method were slightly greater than the field monitoring data, but were much closer to the monitored data than the predictions from other methods. Therefore, the newly proposed method is both applicable to real-world projects and more accurate than existing methods.
Supplemental Material
sj-docx-1-sci-10.1177_00368504221098886 - Supplemental material for Classification of the loosening zones and estimation of the loosening pressures of tunnels in layered jointed rock strata
Supplemental material, sj-docx-1-sci-10.1177_00368504221098886 for Classification of the loosening zones and estimation of the loosening pressures of tunnels in layered jointed rock strata by Zhiwei Wang and Weibin Ma in Science Progress
Footnotes
Acknowledgements
The authors want to express sincere gratitude towards the anonymous reviewers for their through and deep remarks, which have greatly helped improve the paper technically as well as editorially.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was financially supported by the Technology Research and Development Program of the China Railway Corporation (K2021T004), and the Science and technology development plan project of China Academy of Railway Sciences Group Co., Ltd (2020YJ170).
Supplemental material
Supplemental material for this article is available online.
Author biographies
Appendices
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
