An Analytical Model for the Excavation Damage Zone in Tunnel Surrounding Rock

: An accurate theoretical model to predict the extent and mechanical behavior of the excavation damage zone (EDZ) in the surrounding rock of deep-buried tunnel is of great importance for the practical engineering applications. Using the elastic-plastic theory and combined with the analysis on the stress characteristics of the tunnel surrounding rock, this paper present a predict model for the EDZ formation and evolution. A three-zone composite mechanical model was established for the tunnel surround rock and the corresponding stress state and displacement ﬁeld of three zones were derived. The effects of the strain softening and dilatancy during rock deformation was taken into account. The correctness of the proposed model was validated by the existing theoretical models. A sensitivity analysis for different inﬂuencing factors in this model was also performed. The results can beneﬁt for the future numerical and experimental studies.


Introduction
Due to the influence of deep-buried tunnel excavation and high in-situ stress, the hydraulic, thermal, and mechanical properties of the surrounding rock around the excavation face is altered significantly [1], leading to the formation of the excavation damage zone (EDZ) or disturbed rock zone (DRZ) [2][3][4][5][6]. An accurate determination of its extent and mechanical behavior is of great importance for the design and construction of the deep underground engineering.
EDZ has been extensively studied experimentally and numerically. For example, by measuring the ultrasonic wave velocity and acoustic emission [7,8] of the surrounding rock, Falls and Young [9], Meglis et al. [10], and Martino and Chandler et al. [11] determined the formation and development of EDZ during the excavation phase. Even though their real-time performance and high accuracy, experiments cannot be widely applied during the entire design and construction process because they are high-cost and time-consuming. The numerical methods can well complement the experimental approaches. For example, Li and Liu [3] used a two-part Hooke's model implemented into FLAC3D [12] and simulated the EDZ formation and evolution around the ED-B tunnel at the Mont Terri site. The simulation results were found to be highly consistent with those from field tests. Considering the blasting-induced damage, Yang et al. [4] adopted LS-DYNA [13] and simulated the EDZ evolution process along with a deep-buried tunnel excavation by drill and blast. Zhu and Bruhns [14] adopted RFPA2D [15] to model the EDZ of a circular tunnel under a wide range of hydromechanical conditions and the concurrent microdamage evolution was captured. However, due to the complexity of the mechanical behaviour and nonlinear deformation characteristics of the EDZ, it is difficult for the pre-existing numerical models to take into considerations the time dependent dynamic formation and evolution of the EDZ induced by excavation process. And the numerical models are still time-consuming, largely limiting their practical applications. Therefore, an accurate and easy-to-use predicable model is needed to provides convenience for extensive practical engineering design and safety assessment. Some analytical solutions have been proposed for the evaluation of the EDZ development. For example, Schartz and Einstein [16] established an analytical model based on the stiffness ratio between the tunnel support and the surrounding rock. Li and Wang [17] derived an analytical solution for the stress state and deformation field of the supported circular tunnel based on elastic theories and plane strain conditions. Considering the creep effect and the delay installation of the support, Fahimifar et al. [18] proposed the analytical solution to predict the time dependent deformation of the tunnel wall. However, the existing theoretical models are still relatively limited and incapable to comprehensively predict the extent and mechanical behavior of the EDZ, reflecting the omission of some key factors in the mechanical analysis.
To address the current research gap, in this paper, based on elastic-plastic mechanics, a novel analytical model for the EDZ in tunnel surrounding rock was established by further considering the effects of the strain-softening and dilatancy during the surrounding rock deformation. A sensitivity analysis is also performed to investigate the influence of various factors in the model.

Mechanical Division of Tunnel Surrounding Rock
The excavation of the tunnel by drilling and blasting destroys the original stress state of the surrounding rock, and the phenomenon of stress concentration occurs [19][20][21], which leads to the destruction of surrounding rock to form a crushed zone, a plastic softening zone, and an elastic zone. Here, we defined that the range of the EDZ includes the crushed zone and plastic zone. In addition, the lateral pressure coefficient will change as well when the tunnel is under the influence of mining. According to the above analysis, considering the influence factors such as lateral pressure coefficient and support resistance, the elasticplastic mechanical calculation model of the tunnel surrounding rock was established, as shown in Figure 1. formation and evolution of the EDZ induced by excavation process. And the numerical models are still time-consuming, largely limiting their practical applications. Therefore, an accurate and easy-to-use predicable model is needed to provides convenience for extensive practical engineering design and safety assessment. Some analytical solutions have been proposed for the evaluation of the EDZ development. For example, Schartz and Einstein [16] established an analytical model based on the stiffness ratio between the tunnel support and the surrounding rock. Li and Wang [17] derived an analytical solution for the stress state and deformation field of the supported circular tunnel based on elastic theories and plane strain conditions. Considering the creep effect and the delay installation of the support, Fahimifar et al. [18] proposed the analytical solution to predict the time dependent deformation of the tunnel wall. However, the existing theoretical models are still relatively limited and incapable to comprehensively predict the extent and mechanical behavior of the EDZ, reflecting the omission of some key factors in the mechanical analysis.
To address the current research gap, in this paper, based on elastic-plastic mechanics, a novel analytical model for the EDZ in tunnel surrounding rock was established by further considering the effects of the strain-softening and dilatancy during the surrounding rock deformation. A sensitivity analysis is also performed to investigate the influence of various factors in the model.

Mechanical Division of Tunnel Surrounding Rock
The excavation of the tunnel by drilling and blasting destroys the original stress state of the surrounding rock, and the phenomenon of stress concentration occurs [19][20][21], which leads to the destruction of surrounding rock to form a crushed zone, a plastic softening zone, and an elastic zone. Here, we defined that the range of the EDZ includes the crushed zone and plastic zone. In addition, the lateral pressure coefficient will change as well when the tunnel is under the influence of mining. According to the above analysis, considering the influence factors such as lateral pressure coefficient and support resistance, the elastic-plastic mechanical calculation model of the tunnel surrounding rock was established, as shown in Figure 1. As shown in Figure 1, assuming that the surrounding rock is a continuous homogeneous isotropic medium, the original rock stress at the infinite distance of the tunnel is , and the lateral pressure coefficient is . After the tunnel excavation and stress redistribution, the fracture zone, the plastic softening zone and the elastic zone are formed outwardly in the surrounding rock successively. Suppose that the tunnel radius is , the As shown in Figure 1, assuming that the surrounding rock is a continuous homogeneous isotropic medium, the original rock stress at the infinite distance of the tunnel is p 0 , and the lateral pressure coefficient is λ. After the tunnel excavation and stress redistribution, the fracture zone, the plastic softening zone and the elastic zone are formed outwardly in the surrounding rock successively. Suppose that the tunnel radius is R 0 , the external radius of the crushed zone is R b , the external radius of the plastic softening zone is R p , and the support resistance of the tunnel is p i . The elastic theory gives the basic formulas as below [22]:

1.
Yield criterion where σ r is radial stress, σ θ is tangential stress, c is cohesion, and ϕ is internal friction angle.

2.
Balanced differential formula 3. Geometric formula where: ε r is the radial strain and ε θ the tangential strain, and u denotes radial displacement of surrounding rock.

4.
Constitutive formula According to the stress characteristics of elastic zone in the surrounding rock, an equivalent model of elastic zone can be established, as shown in Figure 2. The surrounding rock in the elastic zone can be simplified as follows: the vertical pressure is p 0 , the lateral pressure is λp 0 , the internal pressure is σ ep , and the radial stress at the elastic-plastic interface is σ ep . By using Kirsch's solution [23], the stress distribution in the elastic zone can be obtained as below:

Mechanical Analysis of Three-Zone
where σ r e denotes the radial stress of elastic zone, σ θ e is the tangential stress of elastic zone, R p stands for the radius of plastic zone, and σ ep is the radial stress at the elastic-plastic interface.
When λ = 1, the stress of elastic zone can be calculated through Equation (5) with no consideration for the influence of lateral pressure coefficient.
When r = R p , at the elastic-plastic interface of surrounding rock, there is: At the elastic-plastic interface, the stress also satisfies the Mohr-Coulomb criterion, from which the following can be obtained: Solving Equation (7) gives: Substitute Equation (8) into Equation (6) and the following can be obtained: It can be seen from Equation (9) that, considering the influence of lateral pressure coefficient, the stress state of the elastic zone is not only related to the original rock stress, but also to the lateral pressure coefficient and angle. In other words, the tunnel stress state varies at different angles. The lateral pressure coefficient has a great influence on the distribution of stress field.
where denotes the radial stress of elastic zone, is the tangential stress of elastic zone, stands for the radius of plastic zone, and is the radial stress at the elasticplastic interface. When = 1, the stress of elastic zone can be calculated through Equation (5) with no consideration for the influence of lateral pressure coefficient.
When = , at the elastic-plastic interface of surrounding rock, there is: At the elastic-plastic interface, the stress also satisfies the Mohr-Coulomb criterion, from which the following can be obtained: Solving Equation (7) gives: Substitute Equation (8) into Equation (6) and the following can be obtained:

Displacement Analysis of Elastic Zone
By substituting the stress results into the constitutive Equation (4), the strain of the elastic zone can be obtained as follows: By substituting Equation (10) into geometric Equation (3), the displacement of elastic region can be obtained as follows:

Mechanical Analysis of Plastic Zone
When the external load on the surrounding rock exceeds its strength, the cohesion and internal friction angle of surrounding rock will reduce to varying degrees. This process is known as strain softening. Therefore, the strain softening effect on surrounding rock must be considered in mechanical analysis. For the convenience of calculation, the strain softening model adopted in this paper is shown in Figure 3, in which the cohesive strain softening model is a three line model. Since the internal friction angle changes in a small range and has little influence on the stress distribution, it is assumed that the internal friction angle of rocks in elastic and plastic zones is ϕ 0 , and that in crushed zone is ϕ b . ⎩ 2 ⎭

Mechanical Analysis of Plastic Zone
When the external load on the surrounding rock exceeds its strength, the cohesion and internal friction angle of surrounding rock will reduce to varying degrees. This process is known as strain softening. Therefore, the strain softening effect on surrounding rock must be considered in mechanical analysis. For the convenience of calculation, the strain softening model adopted in this paper is shown in Figure 3, in which the cohesive strain softening model is a three line model. Since the internal friction angle changes in a small range and has little influence on the stress distribution, it is assumed that the internal friction angle of rocks in elastic and plastic zones is , and that in crushed zone is . According to Figure 3, the softening modulus of cohesion could be expressed as: where: denotes the initial cohesion, is the residual cohesion, stands for the critical tangential strain at the interface of elastic zone and plastic softening zone, refers to the critical tangential strain of surrounding rock at the interface between the plastic softening zone and the crushed zone.
where: c 0 denotes the initial cohesion, c b is the residual cohesion, ε θ ep stands for the critical tangential strain at the interface of elastic zone and plastic softening zone, ε θ pb refers to the critical tangential strain of surrounding rock at the interface between the plastic softening zone and the crushed zone.
The surrounding rock of the plastic softening zone is assumed to be incompressible, and ε θ + ε z + ε r = 0. The plane strain condition is ε z = 0 and ε θ + ε r = 0, so the geometric Equation (3) is expressed as follows: To solve the differential Equation (13), we can get u = C 1 r , where C 1 is the integral constant, which is related to the boundary conditions.
From the definition of equivalent strain, the equivalent effect of surrounding rock in the plastic softening zone evolves to be: When r = R p , and ε i = ε θ ep , with the contact conditions of the interface between the elastic zone and the damage zone taken into consideration, the following equation can be obtained: When Equation (15) is substituted into Equation (14), it can be concluded that the effect in the softening zone is as follows: The tangential strain of the plastic softening zone is: The relationship between the cohesion and the radius after the surrounding rock softening can be deduced as follows: At the same time, the dilatancy effect of surrounding rock at yield is considered in the plastic zone as well. The model of dilatancy effect of surrounding rock is shown in Figure 4, where η 1 and η 2 denote the dilatancy gradients of plastic softening zone and crushed zone of surrounding rock respectively, which are assumed to be constant. elastic zone and the damage zone taken into consideration, the following equation can be obtained: When Equation (15) is substituted into Equation (14), it can be concluded that the effect in the softening zone is as follows: The tangential strain of the plastic softening zone is: The relationship between the cohesion and the radius after the surrounding rock softening can be deduced as follows: At the same time, the dilatancy effect of surrounding rock at yield is considered in the plastic zone as well. The model of dilatancy effect of surrounding rock is shown in Figure 4, where and denote the dilatancy gradients of plastic softening zone and crushed zone of surrounding rock respectively, which are assumed to be constant. In the plastic softening zone, the plastic flow rule is as follows: In the crushed zone, the plastic flow rule is as follows: In the plastic softening zone, the plastic flow rule is as follows: In the crushed zone, the plastic flow rule is as follows:

Stress of Plastic Softening Zone
According to Equations (1) and (18), the yield criterion of plastic softening zone is as follows: where: n = 1+sin ϕ 0 1−sin ϕ 0 , m = 2 cos ϕ 0 1−sin ϕ 0 . By substituting Equation (21) into Equation (2), it can be deduced that: To solve Equation (22), we can get: Combined with boundary conditions: r = R p , σ r = σ ep , the following can be obtained: where: 1+n . By substituting Equations (23) and (24) into Equation (21), the distribution of stress field in the plastic softening zone can be obtained:

Displacement of Plastic Softening Zone
According to the flow law and geometric Equation (3) in the plastic softening zone, the following results can be obtained: By solving the differential Equation (26), the following result is obtained Since the boundary condition of the crushed zone is r = R p and u p = u e , it can be deduced that: where: (27), the displacement field of the plastic softening zone can be obtained:

Stress of Fracture Zone
The surrounding rock in the crushed zone not only meets the yield criterion of Equation (1), but also the equilibrium differential Equation (2). The following results can be obtained by combining the two equations: where: With the stress boundary condition being r = R 0 ,σ r b = p i , it can be deduced that: The stress state of the crushed zone can be obtained from Equations (1), (30) and (31):

Displacement of Crushed Zone
According to the plastic flow law and the geometric equation, the following conclusions can be obtained: By solving the differential equation, the result can be obtained With the boundary condition being r = R b and u b = u p ,it can be deduced that: By substituting Equation (35) into Equation (34), the displacement of the crushed zone can be obtained:

Analysis of the Range of Plastic and Crushed Zone
According to the strain softening model, the strength of surrounding rock softens to the residual value at the interface between the plastic zone and the crushed zone. Equation of softening modulus is as follows: On the elastic-plastic interface, it is known from Equation (29): ε θ ep = N. On the interface between plastic softening zone and fracture zone, when r = R b is substituted into Equation (29), the result is as follows: Then: By combining Equations (37)-(39), we can get: At the interface of plastic softening zone and fracture zone, the stress continuity condition are r = R b and σ r p = σ r b . Substituting them into Equations (25) and (32), we can get: By solving the above Equation (41), it can be known that: By combining Equations (40) and (42), the following can be obtained:

Sensitivity Analysis of Influencing Factors
The extent of EDZ directly affects the deformation control and support parameter design of roadway [24]. Equation (43) shows that the extent of EDZ is mainly affected by the original rock stress, lateral pressure coefficient, cohesion, internal friction angle and support resistance. In order to analyze the influence of various factors, the single factor sensitivity analysis method is adopted here. The calculated results are compared with the classical Rubin's solution [25] to verify the applicability of our proposed model.
The Rubin's solution is as follows [25]: Here, the original rock stress (p 0 ) takes the value of 20 MPa, the cohesion (c) is 2 MPa, the internal friction angle (ϕ) is 30 • , the support resistance (p i ) is 0.2 MPa. The calculated results show that the predictions from our model are slightly larger than those of the Rubin's solution. This can be justified by the fact that the factors such as strain softening and dilatancy are further considered in our model. In addition, we use a three-zone composite mechanical model while in the Rubin's solution all the surrounding rock are assumed as ideal elastic-plastic bodies. This indicates that our results are more accurate while Rubin's solution will underestimate the extent of EDZ. The detailed single factor sensitivity analysis is presented as below:

Lateral Pressure Coefficient
When the lateral pressure coefficients are 0.5, 1, 1.5, 2, and 2.5, respectively, and other parameters remain unchanged, the Rubin's solution and the solution in this paper are compared and shown in Figure 5. It can be seen from Figure 5 that the size of plastic zone is directly proportional to the lateral pressure coefficient. The calculated result of plastic zone in this paper is slightly larger than that of the Rubin's solution. When the lateral pressure coefficient exceeds 2, the influence of the lateral pressure coefficient becomes more remarkable. Therefore, the deeper the roadway is buried, the greater the lateral pressure coefficient, and the larger the extent of EDZ. It can be seen from Figure 5 that the size of plastic zone is directly proportional to the lateral pressure coefficient. The calculated result of plastic zone in this paper is slightly larger than that of the Rubin's solution. When the lateral pressure coefficient exceeds 2, the influence of the lateral pressure coefficient becomes more remarkable. Therefore, the deeper the roadway is buried, the greater the lateral pressure coefficient, and the larger the extent of EDZ.

Original Rock Stress
When the original rock stress is set as 5, 10, 15, 20, and 25, respectively, the lateral pressure coefficient is taken as 2, and other parameters remain unchanged, the Rubin's solution of plastic zone distribution extent and the solution in this paper are compared and shown in Figure 6. It can be seen from Figure 5 that the size of plastic zone is directly proportional to the lateral pressure coefficient. The calculated result of plastic zone in this paper is slightly larger than that of the Rubin's solution. When the lateral pressure coefficient exceeds 2, the influence of the lateral pressure coefficient becomes more remarkable. Therefore, the deeper the roadway is buried, the greater the lateral pressure coefficient, and the larger the extent of EDZ.

Original Rock Stress
When the original rock stress is set as 5, 10, 15, 20, and 25, respectively, the lateral pressure coefficient is taken as 2, and other parameters remain unchanged, the Rubin's solution of plastic zone distribution extent and the solution in this paper are compared and shown in Figure 6. It can be seen from Figure 6 that with the increase of the original rock stress, the extent of the plastic zone increases significantly. This is because the shallower surrounding rock releases the stress by fractures and plastic deformation, and the stress concentration zone is transferred to the deeper regions, resulting in the increase of plastic zone range. This can also justify why the superposed mining pressure will increase the plastic zone extent and the difficulty of roadway support. It can be seen from Figure 6 that with the increase of the original rock stress, the extent of the plastic zone increases significantly. This is because the shallower surrounding rock releases the stress by fractures and plastic deformation, and the stress concentration zone is transferred to the deeper regions, resulting in the increase of plastic zone range. This can also justify why the superposed mining pressure will increase the plastic zone extent and the difficulty of roadway support.

Cohesion
When the cohesion is taken as 2, 3, 4, 5, and 6, respectively, and the lateral pressure coefficient is taken as 2, and other parameters remain unchanged, the Rubin's solution of plastic zone distribution range and the solution in this paper are compared and shown in Figure 7.

Cohesion
When the cohesion is taken as 2, 3, 4, 5, and 6, respectively, and the lateral pressure coefficient is taken as 2, and other parameters remain unchanged, the Rubin's solution of plastic zone distribution range and the solution in this paper are compared and shown in Figure 7. It can be seen from Figure 7 that the range of plastic zone decreases significantly with the increase of cohesion. This is because the increase of cohesion leads to the increased strength of surrounding rock, and improved ability of rock mass to resist deformation. In It can be seen from Figure 7 that the range of plastic zone decreases significantly with the increase of cohesion. This is because the increase of cohesion leads to the increased strength of surrounding rock, and improved ability of rock mass to resist deformation. In this case the range of plastic zone can be obviously contained. The integrity of rock mass can be improved by means of grouting or injection anchor, so as to increase the cohesion of rock mass and reduce the extent of plastic zone.

Support Resistance
When the support resistance is set as 0.2, 0.6, 1.0, 1.4, and 1.8, respectively, the lateral pressure coefficient is taken as 9, and other parameters remain unchanged, the Rubin's solution of plastic zone distribution range and the solution in this paper are compared and shown in Figure 8. It can be seen from Figure 7 that the range of plastic zone decreases significantly with the increase of cohesion. This is because the increase of cohesion leads to the increased strength of surrounding rock, and improved ability of rock mass to resist deformation. In this case the range of plastic zone can be obviously contained. The integrity of rock mass can be improved by means of grouting or injection anchor, so as to increase the cohesion of rock mass and reduce the extent of plastic zone.

Support Resistance
When the support resistance is set as 0.2, 0.6, 1.0, 1.4, and 1.8, respectively, the lateral pressure coefficient is taken as 9, and other parameters remain unchanged, the Rubin's solution of plastic zone distribution range and the solution in this paper are compared and shown in Figure 8. It can be seen from Figure 8 that the support resistance of roadway also plays an important role in reducing the extent of the surrounding rock plastic zone. However, the support resistance provided by the support methods such as anchor bolt is very limited, It can be seen from Figure 8 that the support resistance of roadway also plays an important role in reducing the extent of the surrounding rock plastic zone. However, the support resistance provided by the support methods such as anchor bolt is very limited, generally not more than 0.3 MPa. Therefore, for the large deformation roadway, unless the bolt support strength is improved noticeably, the control effect of roadway deformation would be rather poor. Accordingly, the comprehensive support method such as bolt plus grouting should be adopted.

Case Study
Our established theoretical model is used to analyzed the EDZ development of a roadway in a coal mine locating in Pingdingshan, Henan Province, China, subjected to the influence of mining pressures. The extent of EDZ before and after mining was calculated. Since the original cross section of this roadway was a straight-wall semicircular arch, the radius of roadway was firstly corrected by the circumscribed circle radius method: where r * denotes the equivalent radius of roadway; h is the height of the straight wall; and B stands for the net width. In this case, h = 1.4 m and B = 4.2 m, give r * = 2.38 m. The Table 1 below gives the mechanical parameters of the surrounding rock. After mining, the stope line was 100 m away from the roadway, and the roadway was subjected to the supporting pressure of the working face, which is about 5 MPa. In addition, the original stress of surrounding rock increased from 15 MPa to 20 MPa after mining. These parameters were substituted into the Equations (29), (36), (42), and (43), and the results were shown in Table 2, which shows that after mining, the depth of crushed zone increases from 1.1 m to 1.39 m, and the depth of plastic zone increases from 1.73 m to 2.08 m. Therefore, the extent of EDZ further increases after mining, which will significantly shorten the effective anchoring section of the anchor bolt and reduced its anchoring force, and increase the deformation of surrounding rock.

Conclusions
Based on the elastic-plastic mechanics, an analytical model for the excavation damage zone in the tunnel surrounding rock is proposed, which considers the combined effects of the strain softening and dilatancy during the rock deformation. The tunnel surrounding rock is divided into three zones, including the elastic zone, plastic zone, and crushed zone. The corresponding analytical expressions of stress state, displacement field, and extant of EDZ in the tunnel surrounding rock were derived. The comparison and analysis with classic Rubin's solution verify the correctness of the theoretical model. By performing a single factor sensitivity analysis, the influence of different factors such as lateral pressure coefficient, in-situ stress, cohesion, and support resistance was analyzed. The proposed model was further used to evaluate the EDZ development of a roadway subjected to the mining influence. The calculation results show that the maximum displacement of the roadway surface increases from 146 mm to 483 mm, and the depth of plastic zone increases from 1.73 m to 2.08 m The proposed model can not only benefit for the future numerical and experimental studies, but also the extensive engineering practices.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The manuscript data used to support the findings of this study are available from the corresponding author upon request.