Study on the Damage Process and Numerical Simulation of Tunnel Excavation in Water-Rich Soft Rock

: The weakening effect is one of the most important causes triggering large deformation and failure of soft-rock engineering; however, few studies paid attention to damage evolution and constitutive relationship of rock in tensile damage in the excavation unloading and water-weakening process, not to mention the coupling process of unloading and water-weakening. In this paper, the mechanism and engineering characteristics of unloading softening and water-softening of water-rich soft rock are analyzed and summarized. Then with the aid of the strain equivalent principle, the damage of surrounding rock caused by unloading softening and water-softening is coupled, and the compression shear damage and the tensile damage of surrounding rock under the unloading process are analyzed. A damage constitutive model of rock subjected to excavation unloading and water-weakening is proposed considering the inﬂuence of water immersion time, and the proposed model is applied in a newly established ﬁnite element simulation method, which is suitable for excavation in the water-rich soft rock. Based on the mechanical-hydraulic-damage coupled method, the progressive failure process of surrounding rock under the dual softening effects can be reﬂected by the deteriorated parameters of damage elements. Finally, the ﬁeld monitoring data of a typical section in the Xujiadi tunnel is used to verify the applicability and accuracy of the proposed dual softening model and simulation method.


Introduction
Mechanical properties of surrounding rock are not only dependent on its internal factors, such as mineral components, microstructures, and connection types between rock particles, but also significantly affected by the external environment, especially excavation and groundwater [1]. Excavation and water can influence the deformation and strength properties of the surrounding rock so that important differences in mechanical behavior can be expected for surrounding rock under different excavation methods and underground water environments [2]. At present, the weakening effects of excavation and water on the soft rock are widely recognized as important factors inducing engineering disasters in the soft rock strata [3,4]. Consequently, a better understanding of the weakening effect on rock is of urgent importance to the stability and durability of rock engineering.
The surrounding rock softening zone is formed during tunnel excavation and unloading, and the deterioration of the softening zone threatens the safety of the project. Huang et al. [5] studied the deformation, parameters, and fracture characteristics of granite under different unloading conditions and paths and drew the typical stress-strain curves of granite under unloading conditions. The intact limestone and limestone with joints and fissures under different unloading paths, unloading speeds, and initial confining pressures are studied, and the deformation and failure characteristics and parameter deterioration of limestone under various unloading conditions are obtained. It is found that the unloading of rock samples will produce an obvious dilatancy phenomenon, and there are two characteristics of shear failure and splitting tensile failure at the same time [6][7][8][9][10][11][12]. Meng et al. [13] studied the mechanical properties of soft rock under different loading and unloading stress paths by a laboratory test and gave the stress-strain curve of soft rock under loading and unloading and considered that the rock mass is more seriously damaged under unloading condition, and its parameter deterioration degree is higher than that of loading test. Wang et al. [14] discussed the differences of mechanical properties of deep-buried soft rock under different stress paths and found the peak axial strain, ultimate strength and residual strength of soft rock are positively correlated with unloading stress level and unloading rate, and the Mohr-Coulomb strength criterion can better describe the strength characteristics of soft rock than the Hoek-Brown strength criterion. Swansson et al. [15] and Shimamoto [16] both proposed different test methods of the confining pressure unloading test to study the rock strength and fracture development under different confining pressures and unloading paths. Viktorov et al. [17] found that when rock is unloaded after explosion, a considerable decrease occurred in velocity of longitudinal waves in the surrounding rock. Lau et al. [18] developed three laboratory test techniques to quantify the damage by measuring the degradation of elastic properties through periodic unloading-reloading cycles. Amundsen et al. [19] showed that hard rock displays strong tensile fracture characteristics during unloading and inside the rock tensile cracks will be produced. However, most attention is paid to the compression-shear damage; few papers have investigated the model of tensile damage. Furthermore, the Mohr-Coulomb criterion is suitable for compression-shear damage and ignores the tensile damage.
For rock engineering in the water-rich soft rock, except for unloading excavation, underground water immersion is inevitable. A considerable number of researches have been carried out during the past dozen years, to investigate the water weakening effect on mechanical properties of different types of rock. According to the former researches, it is known that the physical-chemical interactions between rock particles and water are greatly influenced by mineralogical composition, grain size, pore volume, and other microstructural properties, which vary significantly among different rock types [20][21][22][23][24][25][26][27][28]. When rock is subjected to external forces or environment changes, microcracks will gradually initiate and propagate and macrocracks will occur, finally inducing rock failure [29][30][31]. The above process, which results in the reduction of rock mechanical parameters from the macroscopical perspective, could also be considered as a damage accumulation process. On one hand, damage models of rock under complex loading or environmental condition have attracted many researchers' attention [30][31][32][33][34]. Nevertheless, few studies focus on the damage evolution during the water-weakening process considering the immersion time. On the other hand, to accurately determine the variation degree of rock mechanical parameters in the process of water-rock interaction, a damage model combining two damage stages induced by unloading and water-weakening should be established.
As for the application of the damage model, Liu et al. [35] established a numerical analysis model of unloading relaxation of rock mass and applied it to the simulation of slope excavation to dynamically simulate the excavation process of a slope. Some articles studied the relationship between post-peak mechanical properties of rock and strain-softening phenomenon, derived the damage evolution equation with softening parameters as variables, and proposed some methods to consider unloading damage effect in finite element numerical calculation [36][37][38]. Lin Peng [39] studied the progressive water inrush caused by particle loss and seepage failure, established a new evolution model of particle loss rate and permeability, and applied it to fluid-solid coupling calculation. Molladavoodi et al. [40,41] used the proposed damage model to simulate of simple loading conditions such as uniaxial tensile and compressive tests and practical and complicated loading conditions such as tunnel excavation. Cho et al. [42] investigated the effect of strain Appl. Sci. 2021, 11, 8906 3 of 28 rate on the damage and failure of tuff through finite element method, which found the tensile strength increases with the increase of strain rate. Armero et al. [43] embedded the finite element into the discontinuous algorithm, established the tensile strength criterion with the maximum normal stress, and analyzed the rock damage and failure process, which is in good agreement with the test. Manouchehrian et al. [44] used an explicit finite element to simulate the unstable rock failures under different unloading conditions. However, most researches about damage simulation focused on the unloading effect; few are based on the unloading and water-weakening coupling damage models. Not to mention considering the time effect in the softening process.
In this paper, the softening stages and softening mechanism of surrounding rock of tunnel engineering in water-rich weak geology are analyzed considering the effects of excavation unloading and underground water. Based on the principle of equivalent strain, dual-softening coupling models of compression-shear damage softening, and tensile damage softening are derived, and the deterioration equations of rock parameters are also derived considering the time effect of water-softening. The applicable conditions of each softening model are summarized in combination with different engineering conditions. A hydro-mechanical-damage multi-field coupling simulation method is proposed based on the dual softening model, which is applied to tunnel engineering to verify its accuracy and applicability.

Staged Softening Process of Water-Rich Soft Surrounding Rock Based on Tunneling
Underground excavation causes a rapid release of in-situ stress in rock and rapid deformation of the surrounding rock to the free surface. The surrounding rock yields due to compression shear or tension, resulting in deterioration of the mechanical properties of the surrounding rock, which is the first softening of the surrounding rock, also called unloading softening. In the water-rich section of geology, unloading softening causes the development of fracture penetration in the surrounding rock, increased permeability, which intensifies the physical, chemical, and mechanical effects of groundwater on the surrounding rock, so that the surrounding rock produces the second softening, that is, water softening.

Undisturbed Stage of Surrounding Rock
After a long time of geological tectonic action, the underground rock is in the original stress state. Due to weathering and seismic action, there are tiny fissures in the rock mass, but the fissures are small in scale and not penetrating. Groundwater cannot flow freely in the fissures, and the softening effect on the rock mass is relatively low and almost negligible. Under this stage, the surrounding rock is in the original stress state, neither unloading softening nor water-softening occurs and all physical and mechanical indexes are relatively stable.

Reinforcement Stage of Surrounding Rock
For tunneling in the water-rich soft rock, forepoling such as pipe sheds or small conduit grouting is usually used to reinforce the surrounding rock. It is generally believed that the pipe shed or conduit spreads and transfers the surrounding rock stress caused by tunnel excavation into the front rock, which plays the role of beam transfer load. Grouting is filling and compacting the rock cracks within the reinforcement layer so that the weak and broken rock forms a relatively complete structure. From a macroscopic point of view, the effect of forepoling is reflected in the improvement of mechanical parameters such as elastic modulus E and shear strength c, ϕ of the surrounding weak and fractured rock. Therefore, to simulate the reinforcing effect of forepoling and anchor rods on the weakly fractured surrounding rock, the equivalent method [45,46] is usually used to improve the physical and mechanical parameters of the rock mass within the influence of grouting.

Unloading Softening Phase during Tunnel Excavation
Tunnel excavation generates massive unloading, the direction of which is mainly related to the shape of the cavern. Unidirectional unloading is for round chambers and vaults of straight-walled arched chambers. Bi-directional unloading is for the location of the sidewalls of straight-walled arches, etc. Most of the highway tunnel cross-sections are in the form of three-centered arches or five-centered arches with a supine arch, which can be considered as unloading in the radial direction only. Figure 1 shows a typical stress-strain relationship for different loading history conditions [12]. There is a clear difference in the mechanical properties of the rock mass during loading and unloading. When the rock mass is loaded or not unloaded, the joint surface of the rock mass is in a closed state. While during unloading, the joint surface slides or opens, resulting in the mechanical parameters of the rock mass being mostly worse than those during loading or not unloading. In general, more unloading leads to more obvious deterioration of the mechanical parameters of a rock mass. Therefore, in the tunnel excavation analysis, the mechanical parameters measured by the loading test will underestimate the softening effect of the unloading of surrounding rock.
elastic modulus E and shear strength , of the surrounding weak and fractured rock. Therefore, to simulate the reinforcing effect of forepoling and anchor rods on the weakly fractured surrounding rock, the equivalent method [45,46] is usually used to improve the physical and mechanical parameters of the rock mass within the influence of grouting.

Unloading Softening Phase during Tunnel Excavation
Tunnel excavation generates massive unloading, the direction of which is mainly related to the shape of the cavern. Unidirectional unloading is for round chambers and vaults of straight-walled arched chambers. Bi-directional unloading is for the location of the sidewalls of straight-walled arches, etc. Most of the highway tunnel cross-sections are in the form of three-centered arches or five-centered arches with a supine arch, which can be considered as unloading in the radial direction only. Figure 1 shows a typical stress-strain relationship for different loading history conditions [12]. There is a clear difference in the mechanical properties of the rock mass during loading and unloading. When the rock mass is loaded or not unloaded, the joint surface of the rock mass is in a closed state. While during unloading, the joint surface slides or opens, resulting in the mechanical parameters of the rock mass being mostly worse than those during loading or not unloading. In general, more unloading leads to more obvious deterioration of the mechanical parameters of a rock mass. Therefore, in the tunnel excavation analysis, the mechanical parameters measured by the loading test will underestimate the softening effect of the unloading of surrounding rock. The stress-strain curve can be divided into four segments, as shown in Figure 2. The OA segment is the elastic stage, indicated by a straight line; the AB part is the unloading yield phase, from which the rock is plastically deformed due to the unloading. The stressstrain relation in this state is a curve, which is in line with the elastoplastic constitutive model, and gradually reaches the peak strength with the increase of the unloading amount; the BC section is the stress-decreasing stage after the peak, and the unloading damage of the rock sample is mainly brittle damage, which can be expressed by linear strain-softening; the CD section is the ideal plastic phase, where the rock reaches the residual strength and shows obvious plastic mobility. The stress-strain curve can be divided into four segments, as shown in Figure 2. The OA segment is the elastic stage, indicated by a straight line; the AB part is the unloading yield phase, from which the rock is plastically deformed due to the unloading. The stressstrain relation in this state is a curve, which is in line with the elastoplastic constitutive model, and gradually reaches the peak strength with the increase of the unloading amount; the BC section is the stress-decreasing stage after the peak, and the unloading damage of the rock sample is mainly brittle damage, which can be expressed by linear strain-softening; the CD section is the ideal plastic phase, where the rock reaches the residual strength and shows obvious plastic mobility.
The deformation of rock samples in the direction of unloading is intense, and brittle damage is obvious during the unloading test [5]. When the confining pressure is large, the shear damage is dominant, and when the confining pressure becomes small, the shear damage gradually transits to tensile damage. During unidirectional unloading, the peak strength, deformation modulus E and residual cohesion c of the rock sample decrease, and the residual friction angle ϕ increases, and the deterioration of rock parameters under bidirectional unloading is greater than that in unidirectional unloading. The cohesive force c decreases because the gelation capacity between fractures is weakened due to the opening of the joint during the unloading process of the rock sample; and the friction angle ϕ increases relatively because of the higher roughness of the resulting tension-shear joint surface. In order to better express the degree of deterioration of mechanical parameters during softening of rock samples, the degree of damage of physical and mechanical parameters is defined as the degree of deterioration. The experimental data of limestone Appl. Sci. 2021, 11, 8906 5 of 28 from actual field projects found that the peak strength of limestone rock samples in the unloaded state deteriorated about 13.9-25.6% and the residual strength deteriorated about 22.8-44.3% compared with the loaded state, and the former was taken as a small value and the latter as a large value when the initial confining pressure was small; the former was taken as a large value and the latter as a small value when the initial confining pressure was large. The deterioration of elastic modulus during unloading is about 7-17%. When the rock is damaged along the fracture surface of the joints, the larger the initial confining pressure is and the slower the unloading rate is, then the greater decrease in elastic modulus is, and vice versa. When the rock sample reaches the residual strength, its cohesion is about 20% deterioration compared with the loading state; and the internal friction angle deterioration is about −12-−5% [10,11]. The deformation of rock samples in the direction of unloading is intense, and brittle damage is obvious during the unloading test [5]. When the confining pressure is large, the shear damage is dominant, and when the confining pressure becomes small, the shear damage gradually transits to tensile damage. During unidirectional unloading, the peak strength, deformation modulus E and residual cohesion c of the rock sample decrease, and the residual friction angle ϕ increases, and the deterioration of rock parameters under bidirectional unloading is greater than that in unidirectional unloading. The cohesive force c decreases because the gelation capacity between fractures is weakened due to the opening of the joint during the unloading process of the rock sample; and the friction angle ϕ increases relatively because of the higher roughness of the resulting tension-shear joint surface. In order to better express the degree of deterioration of mechanical parameters during softening of rock samples, the degree of damage of physical and mechanical parameters is defined as the degree of deterioration. The experimental data of limestone from actual field projects found that the peak strength of limestone rock samples in the unloaded state deteriorated about 13.9-25.6% and the residual strength deteriorated about 22.8-44.3% compared with the loaded state, and the former was taken as a small value and the latter as a large value when the initial confining pressure was small; the former was taken as a large value and the latter as a small value when the initial confining pressure was large. The deterioration of elastic modulus during unloading is about 7-17%. When the rock is damaged along the fracture surface of the joints, the larger the initial confining pressure is and the slower the unloading rate is, then the greater decrease in elastic modulus is, and vice versa. When the rock sample reaches the residual strength, its cohesion is about 20% deterioration compared with the loading state; and the internal friction angle deterioration is about −12-−5% [10,11].

Continuous Water-Softening Stage
During the excavation of underground rock mass, the water-rock interaction includes: (1) the composition and mechanical properties of the rock mass change under the

Continuous Water-Softening Stage
During the excavation of underground rock mass, the water-rock interaction includes: (1) the composition and mechanical properties of the rock mass change under the action of water, (2) the fluid properties and chemical composition of the water itself are also changing. Usually, physical action, chemical action, and mechanical action are complementary and inseparable, especially in the soft rock. As the cementation between diagenetic minerals of soft rock is relatively weak, the gap between mineral particles is large, and the area affected by water is large. Similarly, compared with intact rock, rock mass has more channels conducive to water flow, so rock mass is more influenced by water.
The physical interaction between water and rock includes lubrication and softening effect. Lubrication reduces the cohesion between diagenetic minerals and the friction between rock masses, thus, reducing the shear strength of rock mass. The softening effect is to change the physical properties of rock fracture fillings, and the physical properties of rock mineral composition. According to the definition of the Engineering Geology Handbook [47], the property that decreases in strength after saturation with water is called softening. The rock softening coefficient is defined as K R = σ cw σ c , Where, K R is the softening coefficient of rock; σ cw is the saturation compressive strength of rock; and σ c is the dry compressive strength of rock.
The water-rock chemical interactions include ion exchange, dissolution and solubilization, hydration, hydrolysis, and redox [48]. The water-rock mechanical interactions include pore hydrostatic pressure and seepage force. The hydrostatic pressure reduces the strength of rock mass by reducing the effective stress of the rock discontinuities, and seepage force reduce the strength of the rock by producing a hydraulic fracturing effect on the rock.

Deformation and Damage Characteristics of Surrounding Rock in Water-Rich Geology
In practical engineering, the water-softening effect on rock is usually accompanied by previous unloading softening. Moreover, differential water-rock interaction is caused due to the initial variability of fracture development and its permeability. The surrounding rock in section K14 + 885~K14 + 935 of Xujiadi tunnel is dolomitic limestone. The groundwater in this section is highly developed. After excavation, there is a large amount of water flowing out from the vault in a sprinkling or gushing way, as shown in Figure 3. The rock mass before encountering water can maintain a specific shape. However, after water immersion, the rock mass softens and disintegrates into "melon pieces" of debris, as shown in Figure 4. dry compressive strength of rock.
The water-rock chemical interactions include ion exchange, dissolution and solubilization, hydration, hydrolysis, and redox [48]. The water-rock mechanical interactions include pore hydrostatic pressure and seepage force. The hydrostatic pressure reduces the strength of rock mass by reducing the effective stress of the rock discontinuities, and seepage force reduce the strength of the rock by producing a hydraulic fracturing effect on the rock.

Deformation and Damage Characteristics of Surrounding Rock in Water-Rich Geology
In practical engineering, the water-softening effect on rock is usually accompanied by previous unloading softening. Moreover, differential water-rock interaction is caused due to the initial variability of fracture development and its permeability. The surrounding rock in section K14 + 885~K14 + 935 of Xujiadi tunnel is dolomitic limestone. The groundwater in this section is highly developed. After excavation, there is a large amount of water flowing out from the vault in a sprinkling or gushing way, as shown in Figure 3. The rock mass before encountering water can maintain a specific shape. However, after water immersion, the rock mass softens and disintegrates into "melon pieces" of debris, as shown in Figure 4.

Deterioration Law of Mechanical Parameters of Water-Encountered Rock
The typical unloading stress-strain curve of rock during water-rock interaction is illustrated in Figure 5 [49]. 1 3 1 , , σ σ ε denote the axial stress, confining pressure, and axial strain, respectively. From Figure 5, it can be seen that the unloading stress-strain curves of the rocks are mostly the same at different water-rock interaction times and can be roughly divided into compaction phase, elastic phase, yield phase, damage phase, and residual phase. With the increase of the times of water-rock interaction, the stress-strain curves of rock samples have some obvious trends: the compaction phase gradually grows, the slope of the elastic phase gradually decreases, the yield phase gradually becomes apparent, the peak strength gradually decreases, the axial strain corresponding to the peak strength gradually increases, the drop of stress-strain curve slows down gradually, the brittle damage characteristic gradually weaken, and the residual strength also gradually decreases. It concludes that the rock sample tends to become softer under the water-rock interaction.

Deterioration Law of Mechanical Parameters of Water-Encountered Rock
The typical unloading stress-strain curve of rock during water-rock interaction is illustrated in Figure 5 [49]. σ 1 , σ 3 , ε 1 denote the axial stress, confining pressure, and axial strain, respectively. From Figure 5, it can be seen that the unloading stress-strain curves of the rocks are mostly the same at different water-rock interaction times and can be roughly divided into compaction phase, elastic phase, yield phase, damage phase, and residual phase. With the increase of the times of water-rock interaction, the stress-strain curves of rock samples have some obvious trends: the compaction phase gradually grows, the slope of the elastic phase gradually decreases, the yield phase gradually becomes apparent, the peak strength gradually decreases, the axial strain corresponding to the peak strength gradually increases, the drop of stress-strain curve slows down gradually, the brittle damage characteristic gradually weaken, and the residual strength also gradually decreases. It concludes that the rock sample tends to become softer under the waterrock interaction.
of the rocks are mostly the same at different water-rock interaction times and can be roughly divided into compaction phase, elastic phase, yield phase, damage phase, and residual phase. With the increase of the times of water-rock interaction, the stress-strain curves of rock samples have some obvious trends: the compaction phase gradually grows, the slope of the elastic phase gradually decreases, the yield phase gradually becomes apparent, the peak strength gradually decreases, the axial strain corresponding to the peak strength gradually increases, the drop of stress-strain curve slows down gradually, the brittle damage characteristic gradually weaken, and the residual strength also gradually decreases. It concludes that the rock sample tends to become softer under the water-rock interaction. According to the results of laboratory tests such as water-saturation test, wet-dry cycle test, and field observation, the strength, elastic modulus, cohesion, and internal friction angle of saturated limestone specimens are reduced, and the degree of reduction is related to many parameters such as initial confining pressure, water immersion time of rock samples, osmotic pressure, and internal fractures. In the saturated state, the increase of the confining pressure makes the internal microcracks close, which leads to the strengthening of the deformation resistance of the rock, so the larger parameters such as compressive According to the results of laboratory tests such as water-saturation test, wet-dry cycle test, and field observation, the strength, elastic modulus, cohesion, and internal friction angle of saturated limestone specimens are reduced, and the degree of reduction is related to many parameters such as initial confining pressure, water immersion time of rock samples, osmotic pressure, and internal fractures. In the saturated state, the increase of the confining pressure makes the internal microcracks close, which leads to the strengthening of the deformation resistance of the rock, so the larger parameters such as compressive strength and elastic modulus lead to the smaller deterioration degree. When the osmotic pressure increases, the overall stress-strain curve moves downward, the peak strength, residual strength, and elastic modulus all decrease to some extent, and the degree of deterioration increases. The longer the immersion time is, the smaller the compressive strength, tensile strength, shear strength, and elastic modulus of the rock sample are, and the greater the parameter deterioration, which will reach stability in about 90 days. The higher the number of wet-dry cycles, the lower the compressive strength, elastic modulus and other parameters of the rock sample, and the greater the deterioration degree. Moreover, the more developed the internal fissures and the lower the initial strength lead to the greater deterioration of mechanical parameters of the rock sample in the saturated state.

Damage Evolution Equation of Soft Surrounding Rock Considering Dual Softening of Unloading and Water Encountering
Most of the simulations of underground engineering in water-rich geological sections in literature mainly favor the fluid-solid coupling of stress and seepage fields, but the effects of excavation-induced unloading and water-softening on the stress and seepage fields of surrounding rocks are less studied. Damage is the microstructural changes under certain load and environment that cause the deterioration of the properties of solid materials, and such microstructural changes reaching a certain level will lead to the destruction of solid materials, so in general, the destruction of materials can be regarded as the process of damage accumulation. In this paper, the damage variables are calculated based on the damage degree of the macroscopic physical and mechanical parameters of the rock mass, and the state of the surrounding rock due to excavation unloading and water-softening is defined as the damage field. Figure 6 shows the coupling relationship between the stress field, seepage field, and damage field.
Most of the simulations of underground engineering in water-rich geological sections in literature mainly favor the fluid-solid coupling of stress and seepage fields, but the effects of excavation-induced unloading and water-softening on the stress and seepage fields of surrounding rocks are less studied. Damage is the microstructural changes under certain load and environment that cause the deterioration of the properties of solid materials, and such microstructural changes reaching a certain level will lead to the destruction of solid materials, so in general, the destruction of materials can be regarded as the process of damage accumulation. In this paper, the damage variables are calculated based on the damage degree of the macroscopic physical and mechanical parameters of the rock mass, and the state of the surrounding rock due to excavation unloading and water-softening is defined as the damage field. Figure 6 shows the coupling relationship between the stress field, seepage field, and damage field. The functions in Figure 6 are explained as follows: Action 1 indicates the effect of hydrostatic pressure and seepage forces generated by groundwater on the stress field. Action 2 indicates the effect of changes in the stress field on the seepage field during tunnel excavation. Action 3 indicates the yielding damage of the surrounding rock under the combined effect of surrounding rock pressure, infiltration pressure, and ground stress release. Action 4 indicates the deterioration of physical and mechanical parameters such as strength, elastic modulus, cohesion, and friction angle after damage to the surrounding rock, which in turn causes stress redistribution. Action 5 indicates that the permeability coefficient increases after the damage to the surrounding rock, resulting in changes in the distribution of the seepage field. Action 6 indicates that the damage of the surrounding The functions in Figure 6 are explained as follows: Action 1 indicates the effect of hydrostatic pressure and seepage forces generated by groundwater on the stress field. Action 2 indicates the effect of changes in the stress field on the seepage field during tunnel excavation. Action 3 indicates the yielding damage of the surrounding rock under the combined effect of surrounding rock pressure, infiltration pressure, and ground stress release. Action 4 indicates the deterioration of physical and mechanical parameters such as strength, elastic modulus, cohesion, and friction angle after damage to the surrounding rock, which in turn causes stress redistribution. Action 5 indicates that the permeability coefficient increases after the damage to the surrounding rock, resulting in changes in the distribution of the seepage field. Action 6 indicates that the damage of the surrounding rock is directly caused by the seepage action of groundwater, such as the deterioration of the parameters brought by various physical and chemical actions of groundwater on the surrounding rock.

Evolution Equation of Damage Variables Considering Dual Softening Coupling
In the excavation of the water-rich soft-rock tunnel, according to the analysis in Section 2, the surrounding rock environment is in order of the natural undisturbed state, excavation unloading damage soften state, unloading and water-softening state. Assuming that D 1 is the unloading damage factor of the surrounding rock, D 2 is the water-softening factor, D 12 is the damage factor of two-damage coupling. Figure 7 shows the relationship between damage factors and surrounding rock states.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 30 rock is directly caused by the seepage action of groundwater, such as the deterioration of the parameters brought by various physical and chemical actions of groundwater on the surrounding rock.

Evolution Equation of Damage Variables Considering Dual Softening Coupling
In the excavation of the water-rich soft-rock tunnel, according to the analysis in Section 2, the surrounding rock environment is in order of the natural undisturbed state, excavation unloading damage soften state, unloading and water-softening state. Assuming that 1 D is the unloading damage factor of the surrounding rock, 2 D is the water-softening factor, 12 D is the damage factor of two-damage coupling. Figure 7 shows the relationship between damage factors and surrounding rock states. The rock mechanical properties will change under the dual softening of excavation unloading and water and the change in elastic modulus reflects the degree of damage inside the material. Based on Lemaitre's equivalent strain principle, the damage variable equation for the surrounding rock under unloading and water-softening effects is derived [59]: The rock mechanical properties will change under the dual softening of excavation unloading and water and the change in elastic modulus reflects the degree of damage inside the material. Based on Lemaitre's equivalent strain principle, the damage variable equation for the surrounding rock under unloading and water-softening effects is derived [59]: where, σ 0 , E 0 are the original effective stress and elastic modulus of surrounding rock; σ 1 , E 1 are the effective stress and elastic modulus of surrounding rock after unloading softening; σ 2 , E 2 are the effective stress and elastic modulus of surrounding rock after unloading and water-softening. ε 01 , ε 12 , ε 02 are the strain generated by effective stress σ 1 in the original state of rock, effective stress σ 2 in the unloading softening state of rock, effective stress σ 2 in the original state rock. From Equations (4)- (6), we obtain that: From Equations (7) and (8), we have: Combining Equations (9) and (10), it can be obtained that: Kang et al. [60] and Meng [61] also derived similar equations for water-softening and uniaxial loading coupling. Equation (11) is derived based on the elastic modulus. In this paper, unloading and water-softening coupled effects on cohesion and internal friction angle are also calculated according to Equation (11), although there is a certain error the engineering needs are met.
At present, most of the damage studies of surrounding rock focus on compression shear yield damage, and the study of tensile yield damage is insufficient. Therefore, this paper considers both compression shear and tensile damage modes. Some scholars add the maximum tensile strength criterion to the Mohr-Coulomb strength criterion, trying to make material yield in tension by setting the tensile strength [62]. However, in the simulation, the damage caused by tensile strain is usually the result of compressive stress and excavation unloading, so the tensile strain should be used as an indicator to determine the tensile damage of the surrounding rock.

Water-Softening Damage in Tension
The approach of Mazars to extend the one-dimensional damage to three dimensions is adopted, and the equivalent tensile strain ε is used to replace the uniaxial tensile strain as the softening parameter [35,[63][64][65][66]. Concerning Mazars' approach, the expression for the damage threshold and damage variable D t of the tension element are given as: where, σ tr is the residual tensile strength of rock; ε t is the maximum elastic tensile strain; ε tu is the ultimate tensile strain; E 0 is the initial (undamaged) elastic modulus of rock; ε is the equivalent strain, and ε 1 , ε 2 , ε 3 are the three principal strain components. The angel brackets < > is a function, which is defined as follows: Combining Equations (13) and (14), it can be seen that: when the main strain is tensile strain, it is included in the equivalence strain calculation; when the main strain is compressive strain, it is not considered.
Considering the coupling of unloading and water-softening effects, the tensile damage variable D t12 is deduced using Equation (11) as follows: where, D t1 , D t2 , D t12 are the rock tensile damage variables caused by unloading softening, water-softening, and the coupling of unloading and water-softening, respectively; σ tr,1 and σ tr,2 are the residual tensile strength of the rock in a natural state and saturated state, separately; E 0 is the initial elastic modulus of rock; ε tu is the ultimate tensile strain, which can be derived from the test.

Water-Softening Damage in Compression Shear
Referring to the work of Qin [62], Su et al. [66], it is found that the damage variable is linearly related to the equivalent plastic strain in the surrounding rock after yielding under compression shear. In this paper, the equivalent plastic strain ε p is taken as the softening parameter and the compression shear damage variable D c can be expressed as: where, ε p is the equivalent plastic strain, the expression of which is: where, ε p 1 , ε p 2 , ε p 3 are the three principal strain components of plastic strain, respectively. ε r p is the equivalent plastic strain when rock reaches the residual strength, obtained from the test. D m is the damage value when rock reaches the residual stage.
Considering the coupling of compression shear and water-softening effects, the compression shear damage variable D c12 is derived: where, D c1 , D c2 and D c12 are the compression shear damage variable of rock caused by unloading softening, water-softening, and the coupling of unloading and water-softening, respectively; D m1 and D m2 are the maximum damage variable of rock under unloading softening and water-softening, separately; ε p is the equivalent plastic strain of rock, and ε r p is the equivalent plastic strain when rock reaches the residual strength.

Damage Deterioration Equation of Surrounding Rock Parameters Considering Dual Softening Coupling
As the damage variables of tensile damage and compression shear damage are different, the deterioration of the surrounding rock parameters caused by them is also different.
Generally, after the ultimate tensile strain is reached, the damage variable reaches 1, which means the loss of bearing capacity, so the elastic modulus, cohesion, and internal friction angle parameters will become 0. However, in compression shear damage, there is an obvious strain softening and residual stage of rock, and the mechanical parameters of the surrounding rock when reaching the residual stage is subject to a different degree of deterioration, thus, it is necessary to distinguish their deterioration process and magnitude.
In this paper, after coupling the dual softening of unloading and water, the deterioration equations of the mechanical parameters of rock are derived as follows: 1.
Water-encountered damage under tension 2.
Water-encountered damage under compression shear where, E 0 , c 0 , ϕ 0 are the initial elastic modulus, cohesion, and internal friction angle of rock, respectively; E 1 , c 1 , ϕ 1 are the elastic modulus, cohesion, and internal friction angle of rock in residual stage only considering the unloading softening, separately; E 2 , c 2 , ϕ 2 are the elastic modulus, cohesion, and internal friction angle of rock in residual stage only considering the water-softening, respectively; D E12 , D c12 , D ϕ12 are the damage degree of the elastic modulus, cohesion, and friction angle of the surrounding rock considering the two softening effects; ε p is the equivalent plastic strain; ε r p is the equivalent plastic strain when the rock reaches the residual strength.

Damage Deterioration Equation Considering Time Effect and Dual Softening Coupling
The loosening circle of surrounding rock caused by excavation unloading needs some time to converge and the differential water-rock action makes the areas with more concentrated development and better hydraulic conductivity in the fracture networks, preferentially expand into a hydraulic fracture. Therefore, the unloading softening and water-softening damage of surrounding rock have a time effect, so the deterioration process of the surrounding rock parameters also varies with time. After considering the time effect, the equations of damage variable evolution and rock parameters' deterioration are improved as follows: where, D(ε, t) is the damage variable considering both softening parameters and time, and D(ε) is the damage variable only considering softening parameters, see Section 3.1.
are the time effect functions of deterioration of elastic modulus, cohesion, and internal friction angle, respectively. Surrounding rock unloading softening time is short. This paper assumes that the damage variables in this period are linear, while the water-softening time is longer. Xia Dong [53] studied the deformation and strength parameters of saturated rock samples with water immersion time and concluded that each physical and mechanical parameter of the rock is exponentially related to time. In this paper, assuming that the damage variables of the surrounding rock obey exponential variation with water-softening time, then: where, Y denotes the parameters of surrounding rock, including elastic modulus, cohesion, and internal friction angle, etc. Y 0 denotes the initial parameters of the surrounding rock, Y 1 denotes the final parameters after unloading softening, and denotes the final parameters after water-softening. When t = t 1 , Y=Y 1 . When t = t ∞ , Y=Y 2 . It is easy to obtain A = Y 1 − Y 2 , C = Y 2 , so Equation (29) can be transformed into: where, t = 0 indicates the start of unloading softening of surrounding rock. t = t 1 indicates the end of unloading softening, and the start of water-softening. t = t 2 indicates the end of water-softening. Usually, t 1 is much smaller than t 2 .
According to the damage variable the segmental expression of the damage variable D(ε, t) considering time and dual softening can be obtained as follows: where, B is derived from the test, reflecting the rapidity of damage to surrounding rock under water-softening. The larger the B is, the more rapid the damage to surrounding rock under water-softening, and the faster the damage reaches stability.

Evolution Equation of Permeability Coefficient
In the coupling analysis of stress-seepage-damage, the permeability coefficient of the surrounding rock will vary with the degree of softening damage. The permeability coefficient of porous media can increase by 1-5 orders of magnitude after significant damage occurs [63]. To reflect the influence of the damage field on the permeability of rock mass, a jump factor ξ is introduced to calculate the permeability coefficient of the damaged rock mass.
where, k 0 are the initial permeability coefficient of the rock mass, k D is the permeability coefficient after damage, and D is the rock damage variable.

Applicability of Damage Evolution Equation
To better apply to engineering simulation, the appropriate calculation methods and formulas are classified for simulation of tunnel excavation under different surrounding rock conditions.

Consideration of Unloading Softening
This method of damage variables' calculation applies to the tunnel excavation in weak surrounding rocks that are not influenced by groundwater or where the influence of groundwater is not significant. The corresponding equations and change curves of the unloading damage variables are shown in Equations (34) and (35)

Consideration of Unloading and Water-Softening
This method of damage variables' calculation is suitable for tunnel excavation in w ter-rich soft rock without considering time, which means only the final damage results taken into account. The corresponding equations and change curves of the dual soften damage variables are shown in Equations (36) and (37) and Figure 9.
Tensile damage: , and , are the residual tensile strength of natural and saturated rocks, spectively; and are the maximum damage values of rock under unloading s Tensile damage: Compression shear damage: where, D t and D c are the damage variables in tensile damage and compression shear damage, respectively; σ tr is the residual tensile strength of rock; ε t is the maximum elastic tensile strain; ε tu is the ultimate tensile strain; ε is the equivalent strain; E 0 is the initial elastic modulus of rock; ε p is the equivalent plastic strain; ε r p is the equivalent plastic strain when rock reaches residual strength; D m is the damage variable when rock reaches the residual stage.

Consideration of Unloading and Water-Softening
This method of damage variables' calculation is suitable for tunnel excavation in waterrich soft rock without considering time, which means only the final damage results are taken into account. The corresponding equations and change curves of the dual softening damage variables are shown in Equations (36) and (37) and Figure 9.
, and , are the residual tensile strength of natural and saturated rocks, spectively; and are the maximum damage values of rock under unloading tening and water-softening, separately. Tensile damage: Compression shear damage: where, σ tr,1 and σ tr,2 are the residual tensile strength of natural and saturated rocks, respectively; D m1 and D m2 are the maximum damage values of rock under unloading softening and water-softening, separately.

Consideration of Unloading, Water Softening and Time Effects
This kind of calculation is applicable to the tunnel excavation under water-rich soft surrounding rock and considering the damage deterioration with time, and Equation (38) is used in combination with Equations (36) and (37). This method can not only calculate the final damage value after excavation, but also obtain the rock damage values at any moment during the softening period. The time history curve of damage variable considering dual softening and time effect is depicted with the other two conditions in Figure 10.
where, D 1 , D 2 , and D 12 are the damage variables under unloading softening, watersoftening, and the two softening effects. t = 0 denotes the start of unloading softening, t = t 1 corresponds to the end of unloading softening and the start of water softening, and t = t 2 marks the end of water softening. In general, t 2 − t 1 can be taken as 90 days.
( ) where, , , and are the damage variables under unloading softening, water-softening, and the two softening effects. = 0 denotes the start of unloading softening, = corresponds to the end of unloading softening and the start of water softening, and = marks the end of water softening. In general, -can be taken as 90 days.

Numerical Simulation Method Based on Dual Softening Effects
In ABAQUS finite element analysis, the mechanical effect of groundwater on surrounding rock is typically calculated by the fluid-structure coupling method, and the Mohr-Coulomb constitutive is commonly used for the surrounding rock. However, the damage deterioration process of rock mechanical parameters during excavation cannot be realized using such a method. To address the problem, this section goes beyond the traditional fluid-solid coupling method and proposes a computational method that truly reflects the softening and deterioration process of the surrounding rock after excavation in water-rich soft geology. The specific calculation flow is shown in Figure 11 where our improvements are shown in dashed boxes. The mechanical parameters (strength, modulus of elasticity, cohesion, internal friction angle, etc.) of each element in every calculation step should be deteriorated referring to the calculated results of the previous step in the newly improved simulation method.

Numerical Simulation Method Based on Dual Softening Effects
In ABAQUS finite element analysis, the mechanical effect of groundwater on surrounding rock is typically calculated by the fluid-structure coupling method, and the Mohr-Coulomb constitutive is commonly used for the surrounding rock. However, the damage deterioration process of rock mechanical parameters during excavation cannot be realized using such a method. To address the problem, this section goes beyond the traditional fluid-solid coupling method and proposes a computational method that truly reflects the softening and deterioration process of the surrounding rock after excavation in water-rich soft geology. The specific calculation flow is shown in Figure 11 where our improvements are shown in dashed boxes. The mechanical parameters (strength, modulus of elasticity, cohesion, internal friction angle, etc.) of each element in every calculation step should be deteriorated referring to the calculated results of the previous step in the newly improved simulation method.
The judgment process of damage deterioration is shown in Figure 12. The Mohr-Coulomb criterion is still used, and the damage variables are calculated for the corresponding mechanical parameters after the shear yielding of the surrounding rock so that the parameters' deterioration is executed under compression shear damage. Referring to the analysis in Section 3.1, distinct tensile damage often occurs during the excavation unloading, and the tensile strain is not generated by the tensile stress. The Mohr-Coulomb criterion cannot describe the tensile yield of an element, nor can it represent the tensile yield of an element by setting the tensile strength (usually the tensile stress in simulation is minimal), so it is necessary to judge the tensile damage of the non-compression-shear damage elements by tensile strain, and then perform parameter deterioration for the elements under tensile damage.  The judgment process of damage deterioration is shown in Figure 12. The Mohr-Coulomb criterion is still used, and the damage variables are calculated for the corresponding mechanical parameters after the shear yielding of the surrounding rock so that the parameters' deterioration is executed under compression shear damage. Referring to the analysis in Section 3.1, distinct tensile damage often occurs during the excavation unloading, and the tensile strain is not generated by the tensile stress. The Mohr-Coulomb criterion cannot describe the tensile yield of an element, nor can it represent the tensile yield of an element by setting the tensile strength (usually the tensile stress in simulation is minimal), so it is necessary to judge the tensile damage of the non-compression-shear damage elements by tensile strain, and then perform parameter deterioration for the elements under tensile damage.
Steps for the implementation of the dual softening model based underground excavation analysis in ABAQUS are as follows: first, check the "PEEQ" (equivalent plastic strain) and "E" (strain) in the previous results, extract the node number and the corresponding strain value. According to the process in Figure 12 and the damage evolution equation in Section 3, the damage variable of each element is calculated, and transferred to the corresponding field variables. (Theoretically, the plastic development of each node is different, so the damage degree and post-softening parameters are different). Then copy the softened results to resubmit the next step calculation.

Validation of the Simulation Method Based on Dual Softening Effect
To verify the feasibility and accuracy of the simulation method proposed in this paper, a typical water-rich weak surrounding rock section of the Xujiadi tunnel is adapted to use the flow-solid coupling method and the dual softening model based on stressseepage-damage multi-field coupling, respectively. Finally, the simulation results are compared with the field monitoring data for verification.

Project Overview
The Xujiadi tunnel is a separate highway tunnel in which the left tunnel length is Steps for the implementation of the dual softening model based underground excavation analysis in ABAQUS are as follows: first, check the "PEEQ" (equivalent plastic strain) and "E" (strain) in the previous results, extract the node number and the corresponding strain value. According to the process in Figure 12 and the damage evolution equation in Section 3, the damage variable of each element is calculated, and transferred to the corresponding field variables. (Theoretically, the plastic development of each node is different, so the damage degree and post-softening parameters are different). Then copy the softened results to resubmit the next step calculation.

Validation of the Simulation Method Based on Dual Softening Effect
To verify the feasibility and accuracy of the simulation method proposed in this paper, a typical water-rich weak surrounding rock section of the Xujiadi tunnel is adapted to use the flow-solid coupling method and the dual softening model based on stress-seepagedamage multi-field coupling, respectively. Finally, the simulation results are compared with the field monitoring data for verification.

Project Overview
The Xujiadi tunnel is a separate highway tunnel in which the left tunnel length is 3862 m, the maximum depth of burial is about 538 m; the right tunnel length is 3914 m, the maximum depth of burial is about 539 m. The surrounding rock of K14 + 885~K14 + 935 tunnel section is dolomitic limestone, gray-white with joints and fissures development, serious weathering rock mass, poor cementation, and poor self-stabilization. The basic quality grade of the rock mass is IV 3 . Before excavation, the surrounding rock in this section was relatively intact in appearance, but after excavation, there is a large amount of water flowing out of the surrounding rocks and the drainage volume is large, as shown in Figure 13. Moreover, the surging water took away the fissure filler, making cementing force decrease rapidly, which results in a more serious and frequent rock falling phenomenon, and even produces a local collapse, as shown in and Figure 14. Section K14 + 885 of the Xujiadi tunnel is selected as a typical section for analysis, which features a buried depth of 390 m, and the quality grade of surrounding rock is IV 3 and the top overlying loose body is about 20 m thick. The tunnel was excavated in upper and lower steps, with the upper step excavated first and the initial lining applied in time. After the settlement of the vault and the convergence of the foot, the lower step was removed and the inverted arch was applied in time to maintain the overall stability of the support.

Model Building
Considering the convenience of modeling and computational efficiency, our model only simulates the 20 m surface loose body and 100 m buried rock stratum, and 5670 kPa vertical load is applied at the junction of the loose body and surrounding rock stratum to replace the unmodeled 170 m surrounding rock stratum. The model size and geometric arrangement are shown in Figure 15. The left and right boundaries of the model constrain the horizontal displacement, the bottom boundary constrains the vertical displacement, and the top is the free surface. The tunnel cross-section is in the form of a three-centered circle with a supine arch, with a span of 12.5 m and a height of 10.2 m. The distance between the left, right and bottom boundaries and tunnel is greater than six times the size of the tunnel, so the boundary effect can be ignored.

Model Building
Considering the convenience of modeling and computational efficiency, our model only simulates the 20 m surface loose body and 100 m buried rock stratum, and 5670 kPa vertical load is applied at the junction of the loose body and surrounding rock stratum to replace the unmodeled 170 m surrounding rock stratum. The model size and geometric arrangement are shown in Figure 15. The left and right boundaries of the model constrain the horizontal displacement, the bottom boundary constrains the vertical displacement, and the top is the free surface. The tunnel cross-section is in the form of a three-centered circle with a supine arch, with a span of 12.5 m and a height of 10.2 m. The distance between the left, right and bottom boundaries and tunnel is greater than six times the size of the tunnel, so the boundary effect can be ignored.
The steel arch of initial support is modelled utilizing the equivalent elastic modulus. The load release coefficient of the surrounding rock is used to allocate the proportion of the load borne by the initial support and secondary lining. According to the recommended value of load sharing ratio for the two-lane tunnel in the Highway Tunnel Design Specification (JTG3370.1-2018) [67] and the comprehensive consideration of engineering experience, it is determined that surrounding rock, initial support, and secondary lining bear 10, 30, and 60% of the ground load, separately. The lateral pressure coefficient is taken as 0.8 due to the poor integrity of the weak surrounding rock. In combination with the recommended values in the engineering geological survey report and the Highway Tunnel Design Specification, the material mechanical parameters in the simulation are determined as shown in Table 1 The steel arch of initial support is modelled utilizing the equivalent elastic modulus. The load release coefficient of the surrounding rock is used to allocate the proportion of the load borne by the initial support and secondary lining. According to the recommended value of load sharing ratio for the two-lane tunnel in the Highway Tunnel Design Specification (JTG3370.1-2018) [67] and the comprehensive consideration of engineering experience, it is determined that surrounding rock, initial support, and secondary lining bear 10, 30, and 60% of the ground load, separately. The lateral pressure coefficient is taken as 0.8 due to the poor integrity of the weak surrounding rock. In combination with the recommended values in the engineering geological survey report and the Highway Tunnel Design Specification, the material mechanical parameters in the simulation are determined as shown in Table 1. In the simulation, the permeability coefficient D k of surrounding rock after tunnel excavation is adopted according to Equations (32) and (33). Referring to the physical and mechanical parameters of the site limestone statistics in Sections 2.3 and 2.4, the parameters required for calculating the damage field in the simulation are given in Table 2. Due  In the simulation, the permeability coefficient k D of surrounding rock after tunnel excavation is adopted according to Equations (32) and (33). Referring to the physical and mechanical parameters of the site limestone statistics in Sections 2.3 and 2.4, the parameters required for calculating the damage field in the simulation are given in Table 2. Due to the extremely developed fractures, large groundwater influx, and the obvious disintegration of the surrounding rock at the site, the damage deterioration of each parameter is taken as the larger value in the range.

Simulation Results
The results of the numerical model considering excavation unloading and water softening are compared and analyzed with the results of the fluid-solid coupling in Figures 16-18, which can show the differences and connections between the different algorithms. From Figure 16, it can be seen that, after excavation, the results of both models show a settlement zone above the tunnel arch shoulder and an uplift zone below the tunnel arch base. The uplift is also affected by the Mohr-Coulomb criterion used in the model, which is not discussed here. The results after considering the dual softening show a larger settlement area in the upper part of the tunnel, with a maximum vault settlement of 58.3 mm, which is 23% larger than the result without softening, while the uplift zone in the lower part of the arch is smaller. In the contrast, the non-softening result shows a smaller settlement area in the upper part of the tunnel than that of the softening model, with a maximum vault settlement of 47.5 mm, while the uplift zone in the lower part of the tunnel is larger than the result of softening model. As can be seen from Figure 17, the lateral deformation trend of both models is the same after excavation, and the lateral deformation of the surrounding rock is mostly concentrated at the arch waist on both sides of the tunnel. The horizontal convergence of the non-softening model is 34.0 mm, while the horizontal convergence of the model with dual softening is 40.8 mm, which is 20% larger. Figure 18 illustrates that the plastic zone development trend of both models is roughly the same, and the maximum equivalent plastic strain area both appear in the arch waist position of the tunnel, but the maximum equivalent plastic strain and plastic area considering the dual softening are much larger than the non-softening results, which is mainly caused by two factors. First, due to the physicochemical effect of groundwater on the surrounding rock, the shear strength of the surrounding rock is reduced, making its bearing capacity lower; second, due to the presence of groundwater, the effective stress of the surrounding rock is reduced. According to the Mohr-Coulomb criterion, the stress state of the surrounding rock is prone to reach the ultimate equilibrium state and destabilize. In the computational model, the unloading effect and softening of water immersion after excavation make some of the mechanical parameters of the surrounding rock element, such as elastic modulus, internal friction angle, cohesion, etc., have different degrees of deterioration. The specific degree of damage is expressed as the damage factor calculated by Equations (36)- (38). The mechanical parameters of the surrounding rock under softening damage are calculated using damage factors, which show the damage deterioration of rock parameters. The deterioration of the elastic modulus of rock elements causes the deformation to be larger than that of the undamaged result, and the damage deterioration of the shear strength parameters of the surrounding rock also makes the calculated plastic region of the surrounding rock larger, which is consistent with our simulation results.
after excavation make some of the mechanical parameters of the surrounding rock element, such as elastic modulus, internal friction angle, cohesion, etc., have different degrees of deterioration. The specific degree of damage is expressed as the damage factor calculated by Equations (36)- (38). The mechanical parameters of the surrounding rock under softening damage are calculated using damage factors, which show the damage deterioration of rock parameters. The deterioration of the elastic modulus of rock elements causes the deformation to be larger than that of the undamaged result, and the damage deterioration of the shear strength parameters of the surrounding rock also makes the calculated plastic region of the surrounding rock larger, which is consistent with our simulation results.

Comparison of Numerical Simulation and Field Monitoring Results
In-situ monitoring can keep abreast of the force and deformation of the lining and surrounding rock, according to which the safety of the lining and the stability of the surrounding rock can be judged, and hazards can be predicted in time to avoid construction

Comparison of Numerical Simulation and Field Monitoring Results
In-situ monitoring can keep abreast of the force and deformation of the lining and surrounding rock, according to which the safety of the lining and the stability of the surrounding rock can be judged, and hazards can be predicted in time to avoid construction

Comparison of Numerical Simulation and Field Monitoring Results
In-situ monitoring can keep abreast of the force and deformation of the lining and surrounding rock, according to which the safety of the lining and the stability of the surrounding rock can be judged, and hazards can be predicted in time to avoid construction accidents and ensure the construction safety [68]. In this section, the computational results of different models are compared with the field monitoring results to test the accuracy and applicability of the improved method considering the dual softening in simulating the mechanical response of the surrounding rock after excavation in water-rich and soft geology. Figure 19 shows the layout of pressure sensors in surrounding rock.

Tunnel Vault Settlement and Convergence Displacement Analysis
The monitoring data and simulation results of the vault settlement and horizontal convergence of section K14 + 885 of the Xujiadi tunnel are shown in Figures 20 and 21. The time history curves of vault settlement and horizontal convergence are both showing the trend of "rapid rise-gradual stabilization-rise again-final stabilization", which is related to the bench construction on site. The surrounding rock converges rapidly after the unloading softening of the upper and lower bench excavation, which shows a rapid linear rise in the curve, and then under the softening of groundwater, the surrounding rock deforms slowly until it stabilizes, showing a slow rise after a rapid rise. The final vault settlement is 60.7 mm and the horizontal convergence is 42.7 mm.  According to the comparison of different curves in Figures 20 and 21, when the flowsolid coupling without considering the effects of softening and time is used, the calculated vault settlement and horizontal convergence are 47.5 mm and 33.8 mm, which are smaller than the monitoring results by about 21.7% and 20.8%, respectively. When the calculation model only considers the softening of excavation unloading, the calculated vault settlement and horizontal convergence are 52.5 mm and 36.6 mm, which are smaller than the monitoring results by about 13.5% and 14.3%, separately, which are closer to the field monitoring value than the fluid-solid coupling without considering softening effect. When the simulation simultaneously considers the dual damage of excavation unloading and water-softening without time, the final calculated vault settlement and horizontal convergence are 56.9 mm and 39.9 mm, which are smaller than the monitoring results by about 6.3% and 6.6%, respectively. However, if the simulation does not involve the time variation, then the calculated deformation of surrounding rock only occurs at the excavation moment, which means the surrounding rock deformation is only related to the construction steps and has nothing to do with time. Although the model can calculate the final vault settlement and horizontal convergence more accurately, it cannot reflect the variation of surrounding rock displacement with time during the construction, nor can it predict the appropriate support time during the design and construction. So, when the simulation considers excavation unloading, water-softening and time effect at the same time, the calculated vault settlement and horizontal convergence are 58.3 mm and 40.8 mm, which are about 4.0% and 4.4% different from the monitoring results, respectively, and the error is minimum compared with other cases. The calculated results considering time effect are not greatly different from the results without time effect, but the deformation curve of the surrounding rock calculated by the time-effect model has an obvious trend of "rapid rise-slow growth-rise again-final stabilization", which matches well with the development of the site monitoring data and provides a good reference for the site construction. Figure 22 shows the comparison between the field monitoring data and the calculated values of different models at each measuring point of section K14 + 885. At the vault, According to the comparison of different curves in Figures 20 and 21, when the flow-solid coupling without considering the effects of softening and time is used, the calculated vault settlement and horizontal convergence are 47.5 mm and 33.8 mm, which are smaller than the monitoring results by about 21.7% and 20.8%, respectively. When the calculation model only considers the softening of excavation unloading, the calculated vault settlement and horizontal convergence are 52.5 mm and 36.6 mm, which are smaller than the monitoring results by about 13.5% and 14.3%, separately, which are closer to the field monitoring value than the fluid-solid coupling without considering softening effect. When the simulation simultaneously considers the dual damage of excavation unloading and water-softening without time, the final calculated vault settlement and horizontal convergence are 56.9 mm and 39.9 mm, which are smaller than the monitoring results by about 6.3% and 6.6%, respectively. However, if the simulation does not involve the time variation, then the calculated deformation of surrounding rock only occurs at the excavation moment, which means the surrounding rock deformation is only related to the construction steps and has nothing to do with time. Although the model can calculate the final vault settlement and horizontal convergence more accurately, it cannot reflect the variation of surrounding rock displacement with time during the construction, nor can it predict the appropriate support time during the design and construction. So, when the simulation considers excavation unloading, water-softening and time effect at the same time, the calculated vault settlement and horizontal convergence are 58.3 mm and 40.8 mm, which are about 4.0% and 4.4% different from the monitoring results, respectively, and the error is minimum compared with other cases. The calculated results considering time effect are not greatly different from the results without time effect, but the deformation curve of the surrounding rock calculated by the time-effect model has an obvious trend of "rapid rise-slow growth-rise again-final stabilization", which matches well with the development of the site monitoring data and provides a good reference for the site construction. The calculation error is maximum at the foot of the arch and minimum at the waist of the arch. The results obtained from the fluid-solid coupling differ greatly from the field monitoring results, and the errors of the calculated values are −17.7, −20.7, and −11.7%, separately, with the largest error at the arch waist and the smallest error at the arch foot. model, considering unloading softening, water-softening and time effect, the mechanical parameters of surrounding rock units change accordingly with the construction process, environment, and time so that the final calculation results are in better agreement with the field monitored values. In addition, all the calculated results are consistent with the trend of the pressure distribution monitored in the field, i.e., the surrounding rock pressure at the foot of the arch is greater than the pressures at the waist of the arch and the pressure at the vault is the smallest, which also shows the reasonableness of both calculations. At the inverted, the multi-field coupling model calculated results have about 16.7% error with the monitoring results, while the fluid-solid coupling model calculated results have about 40.4% error. As the invert is subjected to much smaller surrounding rock pressure than other monitoring positions, and typically is not the weak position of the tunnel, it is not discussed in this model. Through the analysis of field monitoring data, fluid-solid coupling results and fluidsolid-damage multi-field coupling results, it can be found that the fluid-solid coupling By comparing the calculated results of different models, the results of the fluidsolid coupling without softening are much smaller than that of new-proposed method considering unloading, water-softening, and time and field monitoring results, which fully indicates that because in the fluid-solid coupling, the physical-mechanical parameters of the surrounding rock element (e.g., E, ϕ, c) do not deteriorated with the softening process. When some elements' parameters should be deteriorated due to softening but still keep the initial condition involved in the calculation, then local yielding damage of surrounding rock will not occur. Therefore, the excavation of rock elements still maintains a good integrity and deformation modulus, resulting in a smaller interaction pressure between the surrounding rock and support than the field monitoring data. On the contrary, in the model, considering unloading softening, water-softening and time effect, the mechanical parameters of surrounding rock units change accordingly with the construction process, environment, and time so that the final calculation results are in better agreement with the field monitored values. In addition, all the calculated results are consistent with the trend of the pressure distribution monitored in the field, i.e., the surrounding rock pressure at the foot of the arch is greater than the pressures at the waist of the arch and the pressure at the vault is the smallest, which also shows the reasonableness of both calculations. At the inverted, the multi-field coupling model calculated results have about 16.7% error with the monitoring results, while the fluid-solid coupling model calculated results have about 40.4% error. As the invert is subjected to much smaller surrounding rock pressure than other monitoring positions, and typically is not the weak position of the tunnel, it is not discussed in this model.

Pressure Analysis of Surrounding Rocks of the Tunnel
Through the analysis of field monitoring data, fluid-solid coupling results and fluidsolid-damage multi-field coupling results, it can be found that the fluid-solid coupling often underestimates the softening of surrounding rock in water-rich soft rock, and overestimates the stability of surrounding rock after excavation, so that the deformation and pressure of surrounding rock after excavation are smaller than the actual situation on site. The new-proposed fluid-solid-damage multi-field coupling method considers the damage field under the combined effect of excavation unloading, water-softening, and time variation, which is better suited for the on-site construction steps and damage phenomena, and the final calculated results are closer to the field monitoring results.

Conclusions
This paper, based on the effect of excavation unloading and water-softening, studied the mechanism of rock damage after underground excavation in water-rich soft rock and derived equations for the evolution of damage variables and the deterioration of rock parameters considering time effect. Based on the above equations, a fluid-soliddamage multi-field coupling method considering the dual softening and time variation was proposed and applied in ABAQUS finite element software. Finally, taking a typical section of the Xujiadi tunnel as an example, the results of different models were compared and analyzed with field monitoring data to verify the accuracy and rationality of the newprosed method for water-rich soft rock excavation. The main conclusions are as follows: When the confining pressure varies from large to small, the rock transitions from compression-shear damage to tensile failure. In the water-rock interaction, different rock fracture developments and their permeabilities cause progressive water-rock interaction. Therefore, the tensile damage and the time effect of water-softening during underground excavation must be taken into consideration.
The fluid-solid coupling analysis ignores the softening damage process of surrounding rock under the action of excavation and water. So, the damage evolution equations of surrounding rock under dual softening actions and different failure modes are proposed. The equivalent plastic strain ε p is used as the softening parameter for compression shear damage, and the compression shear damage variable D c and ε p are approximately linear; the equivalent strain ε is used as the softening parameter for tensile damage, and the tensile damage variable D t and ε are nearly negative linear. From the engineering reality, taking into account the difference between excavation unloading and water-rock interaction, the time segmental function is used to reflect the softening effect of rock damage. Usually, the time of excavation unloading is short, and linear function is used; the time of watersoftening is long, an exponential function is used according to the test results. Meanwhile, the calculation methods of softening damage variables corresponding to different working conditions are summarized.
Following the derivation of damage evolution equations, an improved FEM analysis method based on ABAQUS is proposed for underground excavation in water-rich soft rock. The damage of surrounding rock is judged through equivalent strain and equivalent plastic strain based on the results of the previous step. Then the damage degree and deterred parameters of each damage element are obtained. The deteriorated parameters are used to replace previous parameters for the next calculation to reflect the progressive damage process of the surrounding rock under the dual softening. Finally, the simulated results of different models are compared with the field monitoring results, and it is found that the surrounding rock displacement and pressure calculated by fluid-solid coupling are much smaller than the monitoring data. When dual softening and time effect are considered, the calculated results of surrounding rock deformation and pressure are within 6% error with the field monitoring data, and the calculated time history curve is also close to the trend of field monitoring curve, which proves the reasonableness and accuracy of the method. This paper is also subject to some inadequacies. When the new simulation method is implemented, multiple models need to be built and 'inp' files need to be modified. Although satisfactory accuracy can be obtained by multiple iterations, the manual operation is slightly troublesome. Moreover, the calculation volume increases exponentially with the increase of observed points when the time effect is considered in the simulation. In the future, the process of automatic judgment, replacement, and the iteration of damage field can be realized by programming language to improve the calculation efficiency.