Modeling the Seepage-induced Damage in Water-Rich Fault Zones at Tunnel Faces by Double Point Material Point Methods

During tunnel construction, crossing water-rich fault fracture zones leads to water influx hazards at the tunnel face. Excavation exposes the fault, causing disturbances and changes in the hydraulic gradient. Thus, the tunnel face becomes a low hydraulic head surface. The naturally high permeability of the fracture zone provides a pathway for groundwater flow, resulting in water seepage along the fracture zone and subsequent water influx at the tunnel face. To investigate this interaction, a two-phase double-point material point method was employed aimed at developing a coupled fluid-solid numerical model for excavating submarine fault tunnels. The model incorporates solid-liquid phase interaction, including rock mass permeability, water viscosity, and tunnel geometry. It analyzes the changes in the liquid-phase flow velocity, pore water pressure, and particle trajectories in the fractured zone ahead of the working face after the fracture zone is exposed, revealing the dynamic evolution of water influx following excavation. Additionally, this study discusses the impact of the solid-liquid phase permeability coefficients on the range of water influx hazards. The research findings demonstrate that the two-phase double-point material point method effectively captures the seepage and water influx processes, offering valuable insights into the mechanisms of sudden water influx in tunnel engineering.


Introduction
When tunnel construction traverses challenging geological areas such as water-rich fault zones, construction disturbances can easily trigger tunnel face collapse, sudden water influx, and other severe accidents.These incidents may have a significant impact on both the construction safety and the stability of underground structures in tunnels.
Numerous scholars have investigated the impact of groundwater on the stability of tunnel faces and softening mechanisms in water-rich fractured zones.Anagnostou [1] first utilized various models to demonstrate the effect of groundwater seepage on tunnel face stability.Li et al. [2] developed a stressdamage coupling model for fluid flow channels within fractured zones, whereas Zhou [3] conducted studies on hydraulic fracturing in fractured rock masses and the permeation damage mechanism of cemented bodies within fractured zones.However, in the simulation process, the impact of groundwater is typically represented by factors such as drag forces, viscosity, and hydraulic pressures, and few studies have directly simulated the characteristics of water flow after the exposure of fractured zones during tunnel face excavation.The material point method (MPM), combines the advantages of Lagrangian and Eulerian algorithms, providing an efficient simulation approach for large deformation problems.Cheng [4] employed MPM to simulate tunnel face instability experiments conducted by Mair [5] and provided insights into the deformation characteristics of the tunnel surrounding rock.Fernandez [6] conducted a comparative simulation of shallow-buried tunnel face instability modes using the NLA and MPM.These studies employed a two-phase single-point MPM for simulation.Recently, to better depict the fluid-structure coupling problems between rock, soil, and water, a two-phase double-point MPM has been proposed.This approach utilizes two sets of material points to distinctly characterize solid and liquid phases, enabling a more accurate simulation of their respective motion states.This approach was initially proposed by Bandara et al. [7] and was further developed and integrated into the MPM software Anura3D by Francesca et al. [8].Other researchers utilized this method to investigate dike stability [9] and undersea landslide problems [10].
This study aims to construct a numerical model for tunnel excavation, revealing a water-rich fault fracture zone, based on a two-phase double-point MPM combined with the author's proposed fault cement critical state constitutive model.To analyze the interaction mechanism between underground water and the fault zone under excavation disturbance, the seepage failure process and mud inflow process of the fractured zone affected by groundwater after tunnel face instability, and to discuss the impact of changes in rock strength, permeability, cementation, and other factors on seepage failure during this process.

Two phase double point MPM and FCCS constitutive model
The two-phase double-point MPM, based on the porous media theory assumes that a saturated porous medium is formed by the superposition of two independent continuous media: solid-phase skeleton and the liquid phase.It is represented by two sets of Lagrangian material points, solid material points (SMPs) and liquid material points (LMPs), each of which carries the properties of its respective phase.Two sets of material points exist in three possible domains [11]-a dry soil skeleton containing SMPs, a free water domain containing only LMPs, and a saturated soil body containing both SMPs and LMPs, as shown in Figure 1.A detailed understanding of the two-phase double-point MPM is given in the papers by Bandara [7] and Fern [12].

Mass conservation equations
Since the solid and liquid phases are computed using their respective material points, the masses of all material points remain constant during the calculation process.The conservation equations for solid and liquid masses are as follows: 3

𝐷(𝑛 𝑆 𝜌 𝑆 ) 𝐷𝑡
+     ∇ •   = 0 (2) where   , and   represent the volumetric concentration ratios of the liquid and solid phases, respectively(  equivalent to porosity in saturated soils),   , and   are the densities of the liquid and solid phases, respectively,  is the velocity vector, and D(•)/Dt indicates the material time derivative.
We assume negligible spatial variability in the liquid phase density, i.e.∇  is essentially 0. In addition, solid-phase particles are considered incompressible.The volumetric strain rate for weakly compressible liquids can be derived based on the mass-conservation equation. (3)

Momentum conservation equations
Under any circumstances, the dynamic behavior of a continuous porous medium may be solved by combining the momentum balance equations for the solid and liquid phases.In the material point method, the momentum balance is solved at the grid nodes.During the solution process, the interaction force between the solid and liquid phases is expressed in the form of the drag force   .
= ∇ • ( ′ +     ) +   +     ,       = ∇ • (    ) −   +      (4) where   is the acceleration of the solid phase,   is the acceleration of the liquid phase,  ′ is the effective stress tensor,   is the stress tensor of the liquid phase,  is the body force vector。  is the drag force controlled by Darcy's law, which reflects the interaction between the solid-phase skeleton and the flowing liquid.It depends on the permeability of the solid skeleton, viscosity coefficient of the liquid phase, and relative velocity difference between the solid and liquid phases.However, when the liquid-phase velocity increases, the liquid flow is no longer a nearly stagnant laminar flow.Therefore, an extension of Darcy's law is required, and the formula for calculating the drag force is as follows: = /√    3 (6) In the equations,   is the dynamic viscosity of the liquid ， and   is the liquid intrinsic permeability. is the non-Darcy flow coefficient, The formula for its calculation is shown in Equation (6), primarily referencing the Ergun equation [13], A is set to 150, and B is set at a constant value of 1.75.
Therefore, the intrinsic permeability   ，referenced from the Kozeny-Carman formula [14] is updated as a function of   as shown in Equation (7), where   is the grain diameter Their respective accelerations can be obtained by solving the momentum balance equations for the solid and liquid phases at the grid nodes.These accelerations were used to calculate the node velocities and the corresponding velocities and momenta of the material points, updating the displacement and position of the material points.Node velocities can be used to calculate the strain rates of material points and, in conjunction with constitutive relations, solve for the stress changes in the material points.

Constitutive model of cement soil
Considering the complex mechanical characteristics of water-rich fault fracture zones, this study aims to describe the mechanical behavior of fault rock masses during excavation and the exposure process based on the authors' previously proposed fault cement critical state constitutive model (FCCS).
The fault cement critical-state constitutive model introduces the fault cementation degree parameter (  ) into the modified Cambridge model to characterize the type and degree of cementation between rock particles.The cementation degree (  ) characterizes the strengthening effect of the cementitious materials in the fault and the softening phenomenon caused by cementation dissociation.Considering an undisturbed fault rock mass that has not undergone water-induced softening or shear deformation in its natural state, its cementation degree is denoted as the initial cementation degree  0 .As the fault rock mass underwent water seepage or shear deformation, the cementation between rock particles in the fault rock mass skeleton underwent dissociation damage.The actual cementation degree    weakens from the initial cementation degree  0 , and this process is primarily affected by hydraulic weakening and volumetric yield.
Hydraulic weakening refers to the process of groundwater seepage, wherein soluble substances and fine particles within the cemented material are carried away by seepage, causing movement in the rock mass skeleton.This results in changes in the internal pore structure of the fault rock mass, which simultaneously affects its permeability and degree of cementation.During this process, the volume initially occupied by the cemented material within the rock mass is transformed into a pathway for water flow after hydraulic breakdown, creating a new flow channel.Under the process of seepage-induced degradation of cementation, the cementation degree   of the rock mass decreases, and simultaneously, its intrinsic permeability   increases.Based on the Kozeny-Carman formula, both parameters are controlled by the porosity   .The formula for updating   in with respect to porosity   is derived as shown below: The volume yield, characterized by an excess of shear strain affecting the cementation between faults, leads to the dislocation and shear fracture of crystal grains within the cement.This results in the breakdown of the cemented material, causing a reduction in the shear resistance of the fault rock mass and weakening its characteristics.To represent the cementation weakening caused by volume yield in the fault rock mass, a weakening coefficient  is introduced, which is a function of the plastic shear strain    .
The integration of the two weakening mechanisms yields the calculation formula for the mechanical cementation degree    of the fault, which varies with porosity   and plastic shear strain    as follows: (11) Based on the degree of fault cementation, the cohesive and shear dilation strengthening effects provided by the cementation were expressed separately using the fault mechanical cementation degree    , as illustrated in Equations 12 and 13: In these equations,  ′  represents the cohesive strengthening parameter for cementation, a and b are the fault clay cohesion enhancement coefficients reflecting the degree of cementation enhancement,  ′  is the dilation-strengthening parameter for cementation, and c and d are the fault clay dilation enhancement coefficients reflecting the degree of cementation enhancement of the dilation.The yield function f and plastic potential function g for the fault materials are updated, as shown in Equation 14.
Furthermore, the presence of cementitious materials in the faults enhances the shear strength of the fault rock mass.Therefore, the shear modulus of the fault rock mass   may be simplified as the sum of the shear modulus of the fault skeleton   and the increased shear modulus owing to cementation   , where   can be assumed to be the function of the mechanical cementation degree (   ) In the equation,  represents the Poisson's ratio of the fault skeleton, and  2 is a parameter reflecting the linear increase of shear modulus with cementation degree, with units in kPa.
To smoothly capture the mechanical transition from elastic to plastic behavior in the fault rock mass, this study introduces the secondary loading surface ratio parameter R into the yield surface function.
In the equation,   represents the plastic strain vector, and u is a material parameter controlling the plastic deformation within the yield surface, reflecting the development of plastic strain within the yield surface.A smaller value of u corresponds to a larger plastic strain.If ＞0, the material is in the plastic state; if ＜0, the material is in the elastic state.
When a fault is subjected to an effective external stress (such as excavation unloading), the fault skeleton and cementitious materials are loaded simultaneously.Considering the weakening effect of fault cementation on the elastic stiffness    , the expression for the effective stress in the fault material is as follows: ) Integrating the above equations, the following equations are provided for calculating the stress and strain increments: Based on the UMAT subroutine of the MPM simulation software Anura3D, a computational module for the fault cement critical-state constitutive model was developed in Fortran.This module can be embedded in the Anura3D software to simulate numerical experiments on tunnel excavation, thereby revealing water-rich fault fracture zones.

Tunnel face collapse simulation 3.1. Model test
To analyze the significant deformation of the rock mass and water during the excavation of a tunnel exposing a water-rich fault zone, we constructed a model for tunnel excavation revealing a vertically saturated fractured zone.Four groups of different FCCS constitutive parameters were defined to analyze their impact on the displacement of the solid and liquid phases during the tunnel face collapse process.
A 2D model of the rock tunnel was constructed using the assumed plane strain condition.The model possessed a length of 45 m with a 5-m thick fractured zone in the middle.There were 20 m of moderately weathered granite on either side of the fractured zone.The tunnel was excavated from one side of the granite to expose the fracture zone.The tunnel possesses a diameter of 6 m, a depth of 12 m, and a water body depth of 5 m at the ground surface.A schematic of the model is illustrated in Figure .2. Herein, the fractured zone is represented using the proposed FCCS constitutive model, the moderately weathered granite is described using the Mohr-Coulomb failure criterion, and the liquid phase is described using the Newtonian fluid model.The bottom boundary of the model was fully fixed, and the roller boundaries were set as the left and right boundaries of the model, which were simultaneously applicable to both the solid and liquid phases.Finally, after excavation, a vertical fixity was applied to the upper and lower boundaries of the excavated section of the tunnel to simulate the implemented support measures.Considering computational stability and efficiency, the mesh size of the fractured zone and the excavated section of the tunnel were set to 1 m.Each mesh contains 12 LMPs and 12 SMPs.The fracture zone parameters used in the model are listed in Table .1.For this model, the initial parameters for the FCCS were derived from triaxial test fitting on the exposed fractured zone of the Qingdao Jiaozhou Bay Tunnel.To validate the effects of different parameters on the constitutive model, three additional control parameter groups were established.These parameters involved a certain degree of weakening of the skeleton and cement strength of the fractured zone based on the initial parameters.The aim was to analyze the impact of different parameters on the collapse of the tunnel face.The density, porosity, permeability, and other parameters of the moderately weathered granite and the fractured zone are listed in Table 2. To expedite permeation simulation, the permeability coefficient of the fractured zone in the model was artificially increased.During the simulation process, gravity was initially applied to the entire model to achieve the initial stress equilibrium.Subsequently, the moderately weathered granite in the tunnel area was excavated until the fractured zone was exposed.Fixity was applied to the upper and lower parts of the excavated section of the tunnel to simulate the tunnel lining.No additional support was provided to the tunnel face, allowing free deformation and failure.The collapse distance of the solid-phase particles and flow range of the liquid phase were recorded under different parameter conditions for comparison.

Result comparative analysis
The computational results indicated that significant deformation and instability occurred on the tunnel face only when both the skeletal and cementation strengths of the fractured zone were simultaneously weakened after excavation and exposure.In the other scenarios, the tunnel face experienced only a minor bulging deformation.Furthermore, when the permeability of the fractured zone remained unchanged, alterations in both the skeletal and cementation strengths had minimal impacts on the water inflow process on the tunnel face.Under the four sets of FCCS parameters, noticeable water inflow incidents occurred after the exposure of the fractured zone.
The deformation of the solid phase and its influence range are illustrated in Figure .3(a).It is evident that only the skeleton and cementing softened simultaneously.The model experienced significant instability collapse after the tunnel face was exposed, after 30 s, the impact range reached 14 m.As for models with only weakening skeleton strength, only weakening cement strength, and the model with no weakening in the fractured zone, after being excavated and exposed, only minor bulging deformations occurred on the tunnel face, their deformation values are 0.421, 0.176, and 0.041 m, respectively.When subjected to excavation disturbance, the skeleton bore more deformation, and its strength variation had a more pronounced impact on the deformation of the fractured zone.The water-inflow range is illustrated in Figure .3(b).This indicates that when the permeability coefficient of the fractured zone remains unchanged, the variations in rock strength have no significant impact on the final range of influence of the tunnel water influx.Only in the initial exposure period did the tunnel that experienced a noticeable collapse exhibit increased porosity and permeability owing to rock deformation and fragmentation at the tunnel face.Consequently, within the first 5 s of exposure, the collapsing model demonstrated a significantly high water influx quantity and range compared to the other models.However, due to the relatively constant permeability in most regions of the fractured zone, there was no significant difference in the water influx quantity and range among the different models after 30 s of working face exposure, with ranges varying between 18 and 20 m.
For the model with simultaneous softening of the skeleton and cement, a simple analysis was conducted on the mutual impact of the solid and liquid phases during the collapse process of the tunnel face.The variation curve of the solid-liquid phase influence range during the collapse is depicted in Figure . 4(a), and the solid-liquid phase collapse state diagrams at 2.5 s and 30 s after the tunnel face exposure are shown in Figure . 4(b).When the tunnel face is exposed, the solid-phase strength in the fracture zone plays a dominant role.Insufficient rock mass strength leads to instability after the exposure of the tunnel face, resulting in rapid collapse along an internal sliding surface.After 2.5 s, a stable collapsed body was formed with an impact range of approximately 5 m.During this process, the flow of the liquid phase lagged slightly.However, the water continues to flow forward.After 7.5 s of tunnelface exposure, water flowed out from the toe of the collapsed mass.With increasing water inflow, after 12.5 s of tunnel face exposure, the water inflow reached a level that could mobilize the rock mass, carrying several rock particles forward within the collapsed mass.The subsequent movement of the rock masses was mainly controlled by the water inflow.By observing the development curve of the influence range of the solid and liquid phases, it can be observed that they developed almost in parallel.To better analyze the impact of the tunnel face collapse on its permeability and cement degree, 270 material points within the 5-m zone behind the tunnel face fracture zone were selected.The curves for the average effective porosity, permeability, and degree of cement during the collapse process are illustrated in Figure .5 According to the curve results, the rock mass underwent significant deformation when the tunnel face was exposed.Due to the rotation and movement of skeleton particles, there is a redistribution of pores, and the effective porosity in the unstable area increases from the initial 0.4 to 0.48 within the first second of collapse.In this process, according to the Kozeny-Carman formula, the permeability increases from 7.40710 -9 to 2.110 -8 .Therefore, the development speed of the liquid-phase peak surface in this model was significantly higher, as shown in Figure .3. Simultaneously, with the rotational displacement of the fractured rock mass in this area, the interskeleton cementation quickly underwent volumetric yielding and dissociation.This leads to a sharp decrease in the degree of cement.
Following the initial deformation and preceding the formation of the collapsed body, the scattered rock mass collapsed from its natural cemented state and reaccumulated.During this process, the porosity underwent small-scale fluctuations.Following the formation of the collapsed body, the porosity gradually increased until it was dispersed by the water impact, indicating that slow pore redistribution still occurred within the collapsed body.Upon stabilizing at approximately 0.52, during this process, the permeability increases, and the cement degree decreases.Overall, most pore redistribution and cementation weakening occurred during the collapse.Following the formation of a stable collapsed body, the effective porosity underwent minor fluctuations owing to the impact of the water inflow and the movement of the collapsed body.
In general, a two-phase double-point MPM can effectively capture the large deformation process of two-phase materials and can be applied to simulate practical engineering problems.

Conclusion
This study provides a brief introduction to the principles of two-phase double-point MPM and their solid-liquid phase equilibrium equations.It presents a constitutive model for fractured rock masses considering the effects of cementation.and provides an updated formula for the cement degree based on porosity.Subsequently, vertical water-rich fracture zone models were constructed and simulations were conducted using a two-phase double-point MPM coupled with a fault cement critical state constitutive model.The main conclusions are as follows.1、 The combination of the two-phase double-point MPM and FCCS constitutive models can effectively describe the deformation and collapse processes of the tunnel face after the exposure of a water-rich fracture zone, in addition to the subsequent water inflow process.It provides accurate coupling and individual dynamic descriptions of solid and liquid phases.2、 In the FCCS constitutive model, the skeleton strength plays a predominant role, exerting a significant impact on the deformation and instability of the fracture zone rock mass.However, the contribution of the cementing agents should not be overlooked.The instability of a tunnel face depends on various factors, including the parameters of the fracture zone, ground stress conditions, and water pressure conditions.3、 Once the excavation exposes the fault, the deformation and collapse of the rock mass lead to the redistribution of porosity in the fractured zone.During this process, the porosity of the fractured

Figure. 3
Figure. 3 Solid deformation and liquid influence range curve under different FCCS parameter models.

Figure. 4
Figure. 4 Solid-liquid collapse state diagram at different time points.

Figure. 5
Figure. 5 Curve of average effective porosity, permeability, and cement degree Variation in the fracture zone behind the tunnel face.