Numerical investigation of geostress influence on the grouting reinforcement effectiveness of tunnel surrounding rock mass in fault fracture zones

Xingyu Xu ,Zhijun Wu ,Lei Weng ,Zhofei Chu ,Qunsheng Liu,c ,Yun Zhou

a The Key Laboratory of Safety for Geotechnical and Structural Engineering of Hubei Province,School of Civil Engineering,Wuhan University,Wuhan,430072,China

b Wuhan University Shenzhen Research Institute,Shenzhen,518057,China

c State Key Laboratory of Water Resources and Hydropower Engineering Science,Wuhan University,Wuhan,430072,China

d Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan,430071,China

Keywords: Numerical manifold method (NMM)Grouting reinforcement Geostress condition Fault fracture zone Tunnel excavation

ABSTRACT Grouting is a widely used approach to reinforce broken surrounding rock mass during the construction of underground tunnels in fault fracture zones,and its reinforcement effectiveness is highly affected by geostress.In this study,a numerical manifold method (NMM) based simulator has been developed to examine the impact of geostress conditions on grouting reinforcement during tunnel excavation.To develop this simulator,a detection technique for identifying slurry migration channels and an improved fluid-solid coupling(F-S)framework,which considers the influence of fracture properties and geostress states,is developed and incorporated into a zero-thickness cohesive element(ZE)based NMM(Co-NMM)for simulating tunnel excavation.Additionally,to simulate coagulation of injected slurry,a bonding repair algorithm is further proposed based on the ZE model.To verify the accuracy of the proposed simulator,a series of simulations about slurry migration in single fractures and fracture networks are numerically reproduced,and the results align well with analytical and laboratory test results.Furthermore,these numerical results show that neglecting the influence of geostress condition can lead to a serious overestimation of slurry migration range and reinforcement effectiveness.After validations,a series of simulations about tunnel grouting reinforcement and tunnel excavation in fault fracture zones with varying fracture densities under different geostress conditions are conducted.Based on these simulations,the influence of geostress conditions and the optimization of grouting schemes are discussed.

1.Introduction

The geography of western China is characterized by its complexity,with fault fracture zones commonly encountered during tunnel construction.In these fault fracture zones,the surrounding rock mass usually consists of extensively developed fractures,which leads to a low bearing capacity.As a consequence,the broken rock mass in these zones tends to deform severely,which may lead to large deformation hazards and threatens the safety of tunnel staff(Martino and Chandler,2004;Yan et al.,2015)during the tunneling process.To address large tunnel deformations,various supporting methods(Wu et al.,2022)have been proposed.Among these methods,grouting reinforcement has been widely adopted due to its economic and practical advantages in limiting deformations of tunnel surrounding rock mass,and this limiting effect heavily relies on the migration range of injected slurry(Salimian et al.,2017).Furthermore,the grouting migration process in the broken rock mass is a typical fluid-solid (F-S) coupling process and is heavily influenced by the geostress conditions(H?ssler et al.,1992).Therefore,understanding the influence of geostress conditions on the effectiveness of grouting reinforcement is critical to the optimization of the grouting reinforcement scheme for tunnels in fault fracture zones.

To investigate this influence,numerous laboratory (Sui et al.,2015;Zhang et al.,2017a;Liu et al.,2019a;Zou et al.,2020;Zhou et al.,2021;Weng et al.,2022) and field tests (Guo et al.,2019;Li et al.,2021) have been conducted,and analytical theories (El Tani,2012;Fidelibus and Lenti,2012;Zou et al.,2020) or empirical formulas(Rafi and Stille,2014)have also been proposed to explain or conclude these test observations.Although the above studies did contribute a lot to uncovering the relations between geostress conditions and slurry migration processes(Bunger et al.,2013),due to the difficulties in the representation of the complex geometry and the boundary condition of the broken surrounding rock mass,these obtained analytical theories are difficult to be widely applied.As an alternative,because of the advantages in visibility,resource commitment and characterizing complex stratigraphic conditions,the numerical techniques are widely applied in the investigation of the grouting migration process (Xu et al.,2022).

Generally,the numerical simulation techniques used to investigate the effectiveness of grouting reinforcement can be classified into three types: (1) continuum-based,(2) discontinuum-based and (3) hybrid continuum-discontinuum.Among the three types of methods,the continuum-based methods,such as the finite element method (FEM) (Contini et al.,2007) and finite difference method (FDM) (Abdollahi et al.,2019),are the most developed ones.Based on FEM and FDM,the slurry migration process(Wisser et al.,2005;Kim et al.,2009;Fu et al.,2013) and the mechanical behaviors of grouting reinforced rock mass (Coulter and Martin,2006;Yu et al.,2018) are investigated.However,limited by the continuum-based assumption,the migration process of injected slurry among broken rock masses cannot be explicitly represented,which makes it difficult to accurately evaluate the migration range(Vardakos et al.,2007).To address this limitation,the discontinuum-based methods,which treat the rock mass as an assembly of discrete grains(Li et al.,2022),are adopted to simulate the migration process of injected slurry among the rock fracture networks (Zhou et al.,2017).Among them,the discrete element methods(DEMs)(Saeidi et al.,2013;Zhao and Xia,2018;Peng et al.,2021) and discontinuous deformation analysis (DDA) (Xiao et al.,2019;Nie et al.,2020) are widely used to study grouting related problems.Although the slurry migration process among the broken rock mass can be directly represented,simulating the response of reinforced broken rock masses during the tunnel excavation,which is a typical continuous-discontinuous failure process,is still challenging.

Thus,to better deal with the continuous-discontinuous process,the hybrid numerical method,which combines both advantages of continuum-and discontinuum-based methods (Deng and Liu,2020),has been introduced to simulate problems related to underground tunneling and grouting processes.Among these hybrid methods,the numerical manifold method (NMM) is a promising one and has developed rapidly in recent years.By combining the advantages of FEM and DDA(Zheng et al.,2015;Zhang et al.,2017b;Tan et al.,2022),NMM can directly and easily capture the displacement jumps across discontinuities,large movements of isolated rock masses,and the interactions between rock fractures(Hu and Rutqvist,2020,2022).As a result,NMM has been commonly adopted to simulate continuous-discontinuous processes,such as rock failure(Kang et al.,2022;Ning et al.,2022)and tunnel excavation-related problems (Tan and Jiao,2017;Xu et al.,2019).

However,the application of NMM in grouting reinforcement is relatively limited.Wu et al.(2019) developed an F-S framework based on a zero-thickness cohesive element-based NMM (Co-NMM) to investigate the fluid migration process in rock fractures,which is related to the slurry migration process.Based on this F-S framework,the NMM is extended to simulate the slurry migration process in pre-defined fracture networks (Liu et al.,2021).However,in the above studies,the algorithm of fluid migration algorithm is specially designed to simulate the splitting process induced by the fluid with high injection pressure,while fluid migration in broken surrounding rock mass under relatively low injection pressure is not realistically represented.Additionally,the influence of geostress and the coagulation process of injected slurry are neglected.Therefore,the existing NMM-based methods are not able to simulate the slurry migration as well as estimate the effectiveness of reinforcement on the broken surrounding rock mass under different geostress conditions.

Thus,in this study,to investigate the influence of geostress on the effectiveness of grouting reinforcement,an NMM-based method for simulating tunnel grouting reinforcement (NMMTGR) is developed by incorporating a newly proposed grouting reinforcement algorithm into the Co-NMM based tunnel excavation method.In this newly proposed grouting reinforcement algorithm,to determine the migration network for slurry migration,a detection technique is first proposed to search for rock fractures that connect to the grouting boundaries or existing migration networks.Then,to simulate the slurry migration process under different geostress conditions,an improved F-S framework,which considers the influence of the fracture properties,the states of geostress,and the slurry pressure,has been developed.After that,a bonding repairing algorithm,which simulates the coagulation of injected slurry among fracture networks of surrounding rock mass is proposed based on the zero-thickness cohesive element(ZE)model of Co-NMM.After all these improvements,this newly proposed NMM-TGR is verified by numerically reproducing the slurry migration process in a single fracture and a regular fracture network,and these predicted results are compared with analytical solutions,laboratory observations and the results of previous numerical methods.After the validations,a series of simulations about tunnel grouting and tunnel excavation in fault fracture zones with various fracture densities under different geostress conditions are conducted.Based on these numerical results,the influence of geostress conditions on the effectiveness of grouting reinforcement,and the corresponding optimization of tunnel grouting scheme are discussed.

2.Fundamental principles of using Co-NMM for simulating tunnel excavation

The most novel features of Co-NMM for simulating tunnel excavation are the finite cover systems (Shi,1996),the ZE model(Wu et al.,2018a),and the tunnel excavation algorithm (Wu et al.,2018b).Among these features,the finite cover systems have been introduced in detail by many previous studies(Qu et al.,2014;Yang et al.,2020).Thus,in this section,only the ZE model and tunnel excavation algorithm,which are closely relevant to the newly developed algorithm,are described briefly below.

2.1.ZE model (Wu et al.,2018a)

As shown in Fig.1a,in the simulation of tunnel excavation by Co-NMM,the problem domain is regenerated by the Voronoi tessellation technique (Wong et al.,2018) to better represent the geometry of the broken surrounding rock mass during the formation of large deformation hazards.Then,as shown in Fig.1b,the assemblage of random rock blocks is covered by the triangle mathematical meshes to generate the manifold elements(Ma et al.,2010),which is used to capture the stress and strain of rock blocks.After that,as shown in Fig.1c,ZEs are embedded between rock blocks to capture interactions between two adjacent rock blocks as intact rock.As shown in Fig.2a and b,in this ZE model,the interaction between two adjacent rock blocks consists of tensile stress,σt,shear stress,τ,and compressive stress,which are functions of the opening displacements,o,or sliding displacement,s,in the normal and tangential directions,respectively.In addition,based on theoandsof the ZEs,a failure criterion is adopted to distinguish the commonly encountered three failure types of rock.As shown in Fig.2c,the ZEs are divided into three states,namely the intact,damaged,and failed state,and the ZEs in the failed state can further be divided into three modes,namely Mode I,Mode II and Mode I-II failure.

Fig.1.Generation of the ZE model:(a)Random polygon rock grains of numerical models,(b) Mathematic mesh of NMM,and(c)ZE embedded between two adjacent rock grains.

Fig.2.Constitutive model of ZE (a) Tensile behavior,(b) Shear behavior,and (c) Failure criterion.

2.2.Tunnel excavation algorithm

The simulation of tunnel excavation consists of two steps,namely the application of initial geostress condition and softening excavation.To accurately model the response of surrounding rock mass during the tunnel excavation,it is essential to reproduce the initial geostress condition of the strata precisely.As shown in Fig.3a,to apply a pre-defined geostress condition in the strata model,evenly distributed loads are first added to the boundaries of the strata model.Under the compression of the distributed loads,the deformation and the stress waves are induced in the strata model,which leads the entire model to be in a state of stress imbalance.Thus,to reach the state of geostress equilibrium,a high damping ratio and non-reflection viscous boundaries are introduced to expedite the dissipation of the stress waves.

Fig.3.Tunnel excavation algorithm: (a) Step 1,application of geostress,and (b) Step 2,softening excavation.

After the equilibrium of the applied geostress,the free boundary of the strata model is transformed into the fixed boundary to maintain the deformation of the strata model and prevent the dissipation of the applied geostress.Then,to simulate the excavation of the tunnel,a core softening technique,in which Young’s modulus of the manifold elements inner the excavation section(the blue region as shown in Fig.3b)is reduced step by step,is adopted to simulate the unloading effect induced by the excavation of the tunnel face.By using this technique,the spatial effects of the tunnel face as well as the excavation advance rate can be reasonably taken into consideration.Once the modulus of the elements is lower than a pre-defined value(about one-thousandth of the original modulus value),these softened elements are completely removed,and the tunnel excavation is finished.

3.Grouting reinforcement algorithm based on the architecture of Co-NMM

To enhance the capability of Co-NMM in evaluating the effectiveness of grouting reinforcement,a novel grouting reinforcement algorithm is proposed and incorporated into the Co-NMM.This proposed grouting reinforcement algorithm consists of four parts,namely detection of the migration network,simulation of flowing slurry in the migration network,grouting aperture model,and bonding repair algorithm.

3.1.Detection technique of slurry migration network

In this migration algorithm,as shown in Fig.4a,the slurry is injected through the grouting holes distributed around the tunnel face and flows through migration networks consisting of channels detected among fractures.To identify the slurry migration network,a migration network detection technique is proposed.At the first step of the detection process,an initial migration network,which is connected to the grouting holes is detected (as shown in Fig.4a).Then,this detection technique continues to work throughout the simulation for updating the migration networks.As shown in Fig.4b,three latest formed fractures (lines in green) are detected.Fractures F1 and F2 are connected to the existing migration network (lines in gray),and are added to the migration network,while F3 is isolated by an intact ZE (line in black) and is not connected to the existing migration network,which means F3 cannot be added to the migration network(as shown in Fig.4c).Similarly,by repeating the network detection process until no fracture would be connected with the latest updated migration network,the migration network updating process is completed.

Fig.4.Detecting and updating migration network: (a) Initiation of the migration network,(b) Detection of newly generated cracks,and (c) Updated the migration network.

3.2.Simulation of flowing slurry in the migration network

Based on the detected migration network,which consists of migration nodes and channels as shown in Fig.5a,a linear compressibility model is introduced to calculate the slurry pressure at each node of the migration network based on the inflow rate of slurry at each node:

Fig.5.Simulation of slurry migration: (a) Components of the migration network,(b) Parallel plate assumption,and (c) Bingham slurry model.

wherePn+1andPnare the slurry pressures of a migration network node at the current time step and last time step,respectively;Ksis the bulk modulus of the slurry,Qn+1andSn+1are the total inflow rate and the saturation of a migration network node at the current time step;Δtis the size of a time step,Vn+1andVnis the volume of a migration node at the current time step and last time step,respectively.Vn+1andQn+1can be calculated by introducing a smooth parallel plate assumption and Bingham laminar slurry model (as shown in Fig.5b and c) as follows:

whereViis the volume of migration nodei,Lijis the channel length between nodesiandj;τ0and μ are the yield strength and the viscosity of the Bingham slurry,respectively;Jijis the pressure gradient,agis the grouting aperture,which will be introduced in detail at the next section.

In addition,the calculated slurry pressure is assumed to be linearly distributed along the migration channels (as shown in Fig.6a),and this linearly distributed pressure is transformed into two pairs of equivalent concentrated nodal forces,which are added to the manifold elements(as shown in Fig.6b),as follows:

Fig.6.Application of slurry pressures to the rock mass: (a) Distributed pressure,and(b) Equivalent node force.

whereFiandFjare the equivalent nodal force added at the two ends of the migration channel.

3.3.Aperture model for flowing slurry

According to Eqs.(2)and(3),the behaviors of the flowing slurry are significantly influenced by the grouting aperture,which in turn is influenced by the in-stress state of the surrounding rock mass and the slurry pressure (Zhang et al.,2018).Additionally,the response of the grouting aperture is also influenced by the fracture properties.However,previous studies using NMM-based methods to investigate the grouting reinforcement process have generally ignored these factors (Wu et al.,2019;Liu et al.,2021).Thus,to more accurately simulate the grouting reinforcement process under varying stress conditions,an aperture model is proposed,which considers the influence of rock properties,geostress conditions,and slurry pressures.

In this aperture model,the states of a migration channel (or fracture) are divided into three types (as shown in Fig.7),namely the stress-free state,opening state,and closing state,and the corresponding aperture of different states can be calculated as follows:

Fig.7.Migration channels in the (a) Stress-free state,(b) Opening state,and (c) Closing state.

whereaiis the initial aperture under stress-free conditions,aois the incremental of aperture for the opening channel,anis the reduction of aperture for a closing channel.As for the channel in the stressfree state,since the real fracture surface is not smooth,the serration of the fracture surface provides the space for the slurry migration(as shown in Fig.7a).Therefore,to satisfy the assumption of smooth parallel plates,these rough fractures must be adjusted to an equivalent initial aperture as follows:

where σcis the uniaxial compressive strength and is set to be equal to the joint compressive strength (JCS) (Barton et al.,1985),which gives the assumption that weathering effecting is ignored.In Eq(7),the roughness of fracture surfaces,which is the main factor of the initial aperture,is reasonably taken into consideration by introducing the joint roughness coefficient (JRC).

As for channels in the opening state,the two surfaces of the channel are not parallel in most cases,and the aperture increment may vary along the channels as shown in Fig.7b.Thus,to satisfy the assumption of the parallel plates,the incremental of the channel increment in the opening state must be also adjusted to an equivalent increment as follows:

whereoais the average aperture and can be calculated byoa=（oi+oj）/2 andr=oi/oj,in whichoiandojrepresent the opening displacement at the two ends of the channel,as mentioned in Section 2.1.

As for channels in the closing state (as shown in Fig.7c),to reflect the reduction of aperture,which shows a non-linear deformation response to the compressive normal stress,a nonlinear relation (Bandis et al.,1983) is adopted to calculate the reduction of the fracture aperture under compressive stress:

where σe(MPa)is the effective normal stress acting on the fracture faces,kn0(MPa/mm) andvm(mm) are the initial normal stiffness and the maximum allowable closure,respectively.In addition,to reflect the influence of rock mass property and fracture property,theJRCandJCSof rock fractures are integrated into the calculation ofkn0andvmby adopting the following equations (Bandis et al.,1983):

In addition,when the slurry is injected into closed fractures under a certain injection pressure,the slurry pressure () acts on the two faces of fractures,and shares part of the external compressive stress (σn),which is previously borne by the skeleton of fracture faces alone.Under this circumstance,the part of σnborne by the fracture skeleton is defined as the effective compressive stress σe,which directly affectsan,can be calculated as follows:

By taking place σein Eq.(9)with Eq.(12),the correspondingancan be calculated as follows:

By now,the fracture aperture,af,in three different states can be reasonably represented.However,due to the influence of the fracture roughness on the flowing process of slurry,the geometric fracture aperture does not equal the grouting aperture for the flowing slurry in the migration channels,and should be adjusted to the grouting aperture in advance as follows (Olsson and Barton,2001):

With the above-proposed aperture model for flowing slurry,the influences of rock fracture roughness and rock properties on the response of the grouting aperture to the geostress of rock mass and grouting pressure is able to be reasonably considered by this NMMTGR method.

3.4.Bonding repairing algorithm

With the above approaches,the slurry migration range in the broken surrounding rock mass with different geostress conditions can be numerically determined.Based on this,to model the enhancement of the bearing capacity and integrity of the broken surrounding rock after grouting,a bonding repairing algorithm is proposed based on the ZE framework.

As shown in Fig.8,in this bonding reinforcement algorithm,the ZEs (the gray line segments as shown in Fig.8a) are embedded between two faces of each fracture to simulate the bonding effect of slurry coagulation,and the mechanical strength of these newly embedded ZEs for grouting reinforcement (G-ZEs) are given according to the types of injected slurries.With these G-ZEs,the isolated broken surrounding rock masses(as shown in Fig.8b)are bonded together and behaves as intact rock.However,as shown in Fig.8c and d,unlike the original ZEs,whose two forming boundaries fit well with each other,the two forming boundaries of G-ZEs are subjected to a dislocation,which means that the original constitutive model of original ZEs is not suitable for the G-ZEs.For instance,in Fig.8c,a G-ZE is newly embedded between two boundaries of a fracture in the closing state,and the initial offsets of opening and sliding displacements,oiandsi,occur,which can lead to simulation errors.Under this circumstance,if the original algorithm of ZE is applied directly,the G-ZE will be subjected to initial shear stress,which should not exist in newly bonded fractures,due to the initial sliding displacement,si.Thus,to eliminate the initial stress,an equivalent opening displacement and sliding displacement,oeandse,are proposed to eliminate the initial set of opening and sliding displacementoiandsiby the following equations:

Fig.8.Bonding repairing algorithm: (a) Bonded broken rock blocks,(b) Broken rock blocks,(c) Newly generated G-ZE,and (d) Original ZE.

whereoandsare the calculated opening and sliding displacement,respectively.By using these equivalent opening and sliding displacements,the behavior of G-ZEs can be represented by the original constitutive model of ZEs.

Based on all the aforementioned developments,both the slurry migration process and bonding reinforcement process of broken rock masses in the fault fracture zone can be simulated by Co-NMM,and this newly developed grouting reinforcement algorithm is integrated with the NMM-tunnel method.

4.Numerical validation

In this section,this developed NMM-TGR method is verified through several numerical cases.First,the slurry migration in a single fracture is simulated by the newly proposed NMM-TGR method,and the predicted results are compared with the analytical results.Then,a laboratory test is numerically reproduced by this newly proposed NMM-TGR method and the original NMMbased method,and the predicted results are compared with laboratory test observations.After that,two groups of slurry migration simulations in single fractures under different stress conditions and injection pressures are conducted by the previous NMM-based method,which neglects the influence of geostress,and the proposed NMM-TGR method,respectively.Based on the comparisons,the ability and validity of this NMM-TGR are further verified.

4.1.Descriptions of numerical models and boundary conditions

As shown in Fig.9,the detailed geometry information and boundary conditions of the single fracture model for the following simulations are illustrated.As shown in this figure,a fracture with length,l,of 1.1 m is embedded in the numerical rock samples.The slurry is injected from the grouting point at the left end of the fracture,and the right boundary of the sample is assumed to be impermeable.In addition to the geometric information,all the mechanical parameters of the rock sample and the injected slurry are listed in Table 1.As shown in this table,theJRCandJCSof the embedded fracture are set as 10 MPa and 20 MPa,respectively;the bulk modulus of the injected slurry,Kg,is 2 GPa,the viscosity of the injected slurry,μ,is 1×10-3Pa s and the yield strength,τ0is 10 Pa.

Table 1Mechanical parameters of the rock sample with the embedded fracture and injected slurry.

Fig.9.Numerical rock sample with an embedded fracture.

As for the detailed boundary conditions,the injection pressures,P0,and confining stress,σxand σy,for different numerical cases are listed in Table 2.As illustrated in this table,there are 11 cases of different boundary conditions,and they are classified into 3 groups.Group 1 includes only one case,where σxand σyare 0 and 0.3 MPa,respectively,and theP0is 0.3 MPa.Group 2 includes 5 cases with the same σx(0 MPa) andP0(0.3 MPa) but different σy(0.3 MPa,0.4 MPa,0.5 MPa,0.6 MPa,and 0.7 MPa).As for Group 3,5 cases with the same σx(0 MPa)and σy(0.3 MPa)but differentP0(0.3 MPa,0.25 MPa,0.2 MPa,0.15 MPa and 0.1 MPa)are conducted.

Table 2Detailed boundary conditions for different groups of cases.

4.2.Slurry flow in a single fracture

For Case 1,the distribution of slurry pressure and aperture along the fracture predicted by this improved NMM-Grouting and the previous method are illustrated in Fig.10.For this single fracture model,without consideration of aperture variation,the migration range of the injected Bingham slurry can be calculated as the following equation proposed by Deere and Lombardi (1985):

Fig.10.Predicted results: (a) Slurry pressure distributions and (b) Distribution of grouting apertures along the single fracture.

whereL(m)is the migration range in the single fracture,andPgis the slurry pressure at a certain point along the single fracture.Then,by substituting Eqs.(13) and (14) into Eq.(16),the final migration range under a certain injecting pressure in a single fracture,with the consideration of aperture variations,can be predicted by the following implicit analytical integral equation:

whereP0is the injecting pressure at the grouting point,σnis the normal compressive stress,which is assumed to be uniformly distributed along the fracture and equals σyin this case.According to Eq.(17),the analytical pressure distribution curve can be obtained by formula transformation and is drawn(the blue line)in Fig.10a.In addition,based on Eq.(17),the analytical distribution curves of the grouting aperture,ag,along the fracture can be deduced and calculated by the following equation:

This analytical distribution curve of the grouting aperture is also drawn(the blue line) in Fig.10b.

As shown in Fig.10a,for the result predicted by the previous method,the slurry migration range is 1.065 m and the slurry pressure,Pg,decreases linearly with a constant gradient with the increase of migration distance,x.However,for the results predicted by the proposed NMM-grouting,the slurry migration range decreases to 0.881 m and the slurry pressure decreases nonlinearly with the increase of migration distance,x.This difference comes from the different distribution of grouting apertures along the fracture.As shown in Fig.10b,without taking stress conditions and slurry pressure into consideration,the grouting aperture predicted by previous methods keeps constant along the rock fracture.However,by taking stress conditions and slurry pressure into consideration,the grouting aperture decreased dramatically from 1.414×10-4m to 0.975×10-4m under the confining stress,σy.In addition,the predicted distribution of slurry pressure and grouting aperture along the fracture agree well with the analytical solutions,which validates the accuracy and reliability of the simulation results using this proposed NMM-Grouting method.

4.3.Slurry flow in fracture networks

In this section,to further validate the ability of this proposed NMM-grouting method in simulating slurry flow in fracture networks,a laboratory test performed by Modellf?rs?k (1987),which has been widely used to validate the numerical grouting program,is reproduced by the improved NMM-TGR method and the original NMM-based method.As shown in Fig.11a,the laboratory test is conducted in the tessellation,which consists of two plates of plexiglass with dimensions of 1200 × 1000 × 15 mm.In this tessellation,a lattice channel network of 1×5 mm lattice channel is pre-set in a symmetrical pattern between the two plates.The slurry is injected through the grouting point in the center of the tessellation with a constant hydraulic pressure head 0.48 m (around 4704 Pa).

Fig.11.Modelling migration of slurry in a fracture network: (a) Experimental tessellation and (b) Numerical model.

To reproduce this laboratory test numerically,a numerical model containing 120,000 square grains with an edge of 1 mm is established,as shown in Fig.11b,which is discretized by two sets of parallel fractures with an initial aperture,ao=1 mm,and a constant fracture spacing of 0.1 m.Young’s modulusEand Poisson’s ratio of the plexiglass is 3 GPa and 0.3,and theJRCof the channel wall is set as 1.The bulk modulus,viscosity,and yield strength of the grout slurry areKs=2000 Pa,μ=0.035 Pa s and τ0=3 Pa,respectively.The constant injection pressure in the injection point isP0=4707 Pa.As shown in Fig.12,the predicted slurry distribution in this migration network at four time points,namely 2.5 s,6 s,22 s,and 1426 s,are compared with the corresponding experimental results.As shown in this figure,the predicted slurry migration ranges by the two numerical methods agree well with the experimental observations (as shown in Fig.12a),which demonstrates that the proposed NMM-TGR method can successfully simulate the migration process in the fracture networks.In addition,as shown in Fig.12b and c,the predicted fluid migration process by the two methods shows no obvious difference when the grouting process last for a short time,such as 2.5 s,6 s,and 25 s.However,when the grouting process last for a long time,such as 1426 s,an apparent difference between the numerical results predicted by the improved NMM-TGR method and the original NMM-based method can be observed.This shows the numerical result predicted by the improved NMM-TGR method is closer to the laboratory results compared with the original NMM-based method,demonstrating the advantages of the proposed NMM-TGR method on simulating the slurry migration process.This advantage comes from the following reason.As shown in Table 3,based on Eqs.(7),(8) and(15),the predicted grouting apertures are 1 mm and 1.01 mm by the original grouting algorithm and improved grouting algorithm of NMM-TGR,respectively.Thus,due to a more accurate characterization of grouting apertures,the NMM-TGR method can simulate the slurry migration process more accurately.

Table 3Apertures of migration channel in the migration process simulated by the original and improved algorithm of NMM.

Fig.12.Comparison between(a)Experimental observations,(b)Numerical results predicted by the improved NMM-TGR method,and(c)Numerical results predicted by the original NMM-based method.

4.4.Effects of stress on the slurry migration process in a single fracture

The effects of stress on the grouting migration are investigated in Cases 2 to 6 of Group 2.As shown in Fig.13a,all five distribution curves of slurry pressure predicted by the previous method completely coincided with each other,which is led by the incapacity of the previous method for characterizing the influence of stress.However,for the results predicted by the proposed NMM-TGR method,the migration ranges decrease from 0.881 ×10-4m to 0.561 ×10-4m with σyincreasing from 0.3 MPa to 0.7 MPa (as shown in Fig.13a).In addition,as shown in Table 4,with the increase of σy,the errors increased from 21% to 91%,which indicates that neglecting the influence of the stress condition on the grouting aperture can lead to a seriously overestimated migration range.Meanwhile,as shown in Fig.13b,for the cases conducted by the previous method,the grouting aperture along the fracture was kept constant with the initial aperture (1.414 × 10-4m).However,for results predicted by NMM-TGR,the maximum grouting aperture of each case occurs at the grouting point and the minimum grouting aperture of each case occurs at the edge of the migration range.With the increase of σy,the maximum grouting apertures of the five cases decrease from 1.414 × 10-4m to 0.869 × 10-4m and the minimum apertures of the five cases decrease from 0.975×10-4m to 0.639 ×10-4m.

Table 4Slurry migration range under different levels of stress.

Fig.13.Effect of stress on the (a) Slurry pressure distributions and (b) Migration apertures along the single fracture.

4.5.Effects of injection pressure on slurry migration process in a single fracture

The effects of the injection pressure on the grouting migration are investigated by the five cases of Group 3.As shown in Fig.14a,the migration range predicted by the previous methods decreases with the injection pressures linearly.However,for the cases simulated by the NMM-TGR,the migration ranges are smaller and decrease non-linearly along the fracture.In addition,as shown in Table 5,with the decrease of the injection pressure,the errors of predicted migration ranges increase from 21% to 40%,indicating that the neglect of the influence of slurry pressure on the grouting aperture can also lead to an overestimated migration range.But this overestimation induced by neglecting the influence of injected slurry pressure is much smaller than that induced by neglecting the influence of stress conditions.As shown in Fig.14b,the distributions of the grouting aperture along the fracture predicted by the previous methods remain the same,while the maximum apertures of 5 cases predicted by the improved NMM-grouting decrease with the injection pressure,and the grouting aperture gradually decrease to the same initial minimum grouting aperture.

Table 5Slurry migration range under different levels of injection pressures.

Fig.14.Effect of slurry injection pressures on the (a) Distribution of slurry pressure and (b) Distribution of grouting apertures along the single fracture.

5.Influence of geostress on the grouting reinforcement effectiveness

To further validate the capability of this proposed NMM-TGR method in evaluating the effectiveness of grouting reinforcement and investigate the influence of geostress on it,simulations of the slurry migration process and the tunnel excavation process in a fractured deep stratum model are conducted.To achieve the simulations,a stratum model of 60 × 60 m (as shown in Fig.15) is established and discretized into 46,984 broken isolated rock blocks,which are of a mean size of 0.15 m and contain 430,892 elements.In the center of the fault fracture zone,a tunnel (the circle region in the yellow)with a diameter of 6 m is placed in advance to simulate the excavation process.Then,12 grouting points are arranged on a circle with a radius of 4.5 m around the tunnel,and a constant slurry injection pressure of 0.5 MPa is applied to simulate the slurry migration process.On the outer boundary of the model,12 sets of different distributed loads are added in both horizontal (σh) and vertical(σv)directions to reproduce 12 cases with different states ofgeostress.Details about the states of the geostress for the 12 cases are listed in Tables 6 and 7.Finally,4 sets of stratum models with different fracture densities establish to investigate the influence of geostress on the grouting reinforcement effectiveness in fault fracture zone with different fracture densities,and details about fracture densities for the 4 sets of models are listed in Table 13.

Table 66 numerical cases with different levels of geostress.

Table 76 numerical cases with different coefficients of lateral stress.

Fig.15.Model of fault fracture zones.

After the construction of the geometric model and determination of different boundary conditions in different cases,appropriate parameters of surrounding rock mass,injected slurry and the reinforced rock mass are also carefully selected.For parameters of broken surrounding rock mass,a typical soft rock mass,which is a kind of mudstone,in the fault fracture zone in the Bayu tunnel of the Sichuan-Tibet Railway (Xue et al.,2020) is adopted for calibration.Then,a trial-and-error process,where the mechanical parameters are calibrated by the comparison between numerical results and laboratory results (Zhang,2012 (in Chinese)),is conducted to obtain the detailed mechanical parameters(as shown in Table 8) of the fractured rock mass.As shown in Fig.16,with the calibrated mechanical parameters,the typical failure patterns of soft rock mass in uniaxial compression and Brazilian tensile test are well reproduced.In addition to these predicted failure patterns,as shown in Table 9,the predicted numerical results,namely,UCS and tensile strength of BD tests,elastic modules and Passion’s ratio,also fit the laboratory results well,with errors within 4%,which further proves the validity of these mechanical parameters.The parameters of injected slurry,which is a typical cement slurry,were obtained from a reference (Liu et al.,2016) and the corresponding parameters are listed in Table 10.Then,the mechanical parameters of reinforced rock mass are determined based on the laboratory results from a previous reference (Liu et al.,2016),where the parameters of newly generated ZEs are increased by 100% for shear stiffness and 93% for shear peak strength respectively.And the detailed mechanical parameters of the grouting reinforced rock are listed in Table 10.

Table 8Calibrated Mechanical parameters.

Table 9Comparison of predicted macromechanical properties with laboratory results.

Table 10Mechanical and grout slurry parameters of the deep tunnel model.

Fig.16.Predicted laboratory results: (a) Brazilian tensile test and (b) Uniaxial compression test.

Finally,by using these mechanical parameters of broken isolated surrounding rock mass,injected slurry and reinforced rock mass,the excavation process and grouting reinforcement process are simulated in the models with 4 sets of facture densities under 12 sets of geostress conditions,and the predicted results are illustrated and discussed below.

5.1.Influence of geostress levels on the grouting reinforcement effectiveness

Due to the influence of the stratigraphic movement and different burial depths of the strata,different levels of geostress are frequently encountered during tunnel construction.In this section,6 cases with varying levels of geostress are conducted to investigate the corresponding influence on grouting reinforcement effectiveness.The numerical results are illustrated from two perspectives.Firstly,the influences on the migration process of the injected slurry in the broken surrounding rock mass are illustrated in Fig.17 and Table 9.Secondly,the influences on reinforcement effectiveness are illustrated in Fig.18.

Fig.17.Slurry migration range in fractured strata models with different levels of geostress: (a) 0.5 MPa (b) 1 MPa (c) 2 MPa (d) 4 MPa (e) 8 MPa,and (f) 16 MPa.

Fig.18.Comparison of convergence in tunnels without and with grouting reinforcement under different levels of geostress:(a)Case 1,(b)Case 2,(c)Case 3,(d)Case 4,(e)Case 5,and (f) Case 6.

As shown in Fig.17,the NMM-TGR is capable of distinguishing the different migration ranges in the fault fracture zone with varying levels of geostress.As shown in Fig.17a and b,it can be observed under low geostress conditions,the total area of the slurry migration range (Amr) is larger,which implies that a larger reinforcement range and a better reinforcement effect can be achieved.However,due to the low level of geostress,the migration range of a single grouting hole (Amrs) is also quite large,resulting in a significant overlap of migration range (Ao) and a low effective grouting rateREG,which is defined as follows:

whereAsmrsis the sum ofAmrs,andNis the number of grouting holes.As listed in Table 11,with the increasing level of geostress from Case 1 to 4,theAmrsof each grouting hole decreases from 78.54 m2to 7.55 m2,and theAmrof each case decreases from 262.52 m2to 69.95 m2.However,with the increasing level of geostress,Aodecreases from 679.96 m2to 20.74 m2,andReincreases from 27.86% to 77.13%.When the level of the geostress increases to 8 MPa and 16 MPa,Amrsfurther decreases to 2.62 m2and 1.36 m2,and gaps begin to appear between the grouting migration islands of adjacent grouting holes.This indicates that a complete and continuous grouting reinforced zone can no longer be achieved when the value of geostress is larger than 8 MPa (as shown in Fig.17e).However,in the cases with high levels of geostress,very highRe,almost reaching 100%,can be achieved,and very lowAoof 0 can also be achieved.

Table 11Slurry migration ranges of numerical Cases 1 to 6.

For the predicted grouting reinforcement effectiveness,Fig.18 shows that the differences in tunnel wall convergence under different levels of geostress are well captured by this improved NMM-TGR method.As shown in Fig.18a-c,for the bare tunnels without grouting reinforcement,the convergence of the tunnel walls accelerates with increasing levels of geostress.At step 700,000,the surrounding rock mass deforms to varying degrees,and the section area of the tunnels is reduced from 28.27 m2to 20.67 m2,19.52 m2and 18.43 m2,respectively.However,for the cases with grouting reinforcement,due to the grouting reinforced zone (the zone in gray) around the tunnel,the deformation of the surrounding rock mass is limited and can even be ignored.With the further increase in the geostress level(4 MPa),as shown in Fig.18d,the grouting reinforced zone fails to cover any of the surrounding rock mass around the tunnel,leaving a layer of broken isolated rock mass on the inner wall of the tunnel.As a result,at step 700,000,the broken isolated rock on the top and the upper left part of the inner wall is spalling into the tunnel.However,even though the unreinforced rock of the inner wall is broken and spalled,the grouting reinforced zone remains intact and inhibits the initiation of fractures.Thus,the deformation of the tunnel’s inner wall is still well-controlled in Case 4,where the area of the tunnel section is 15.23 m2and 22.94 m2in the tunnels without and with grouting reinforcement,respectively.For Case 5,due to the higher level of geostress (8 MPa),the grouting reinforced zone is no longer a complete ring but consists of several isolated reinforced islands with large gaps in-between,leading to a sharp decrease in grouting reinforcement effectiveness.Additionally,the higher level of geostress leads to the initiation of fractures in the isolated grouting reinforced islands,further weakening the grouting reinforcement effectiveness.Due to the influence of the increasing level of geostress,the difference in the area of the tunnel section between the cases without and with grouting reinforcement decreases to 5.12 m2,which is 7.71 m2for Case 4.For Case 6 with a higher level of geostress (16 MPa),it is evident that the range of the grouting reinforced zone further decreases,and the isolated reinforced islands are separated by larger ranges of the broken surrounding rock mass.Moreover,these isolated reinforced islands fracture into broken rock masses again during the excavation process.Thus,as shown in Fig.18f,at step 700,000,the area of the deformed section of the tunnels with and without grouting reinforcement is 4.64 m2and 4.21 m2,respectively,indicating that the effectiveness of grouting reinforcement is almost lost.

5.2.Influence of lateral geostress ratio on the grouting reinforcement effectiveness

Besides the level of geostress,the effectiveness of grouting reinforcement on the broken rock mass in fault fracture zones is significantly affected by different lateral geostress ratios.Therefore,another 6 cases are simulated to investigate the corresponding influence of different lateral geostress ratios on the grouting reinforcement effectiveness.

For the migration process of these 6 cases with different lateral geostress ratios,it is evident that the migration range changes not only in size but also in shape as σhchanges.For instance,in Case 7,where the lateral geostress is 2 MPa,the migration range of a single grouting hole is in an ellipse shape with a long axis of 2.1 m in the vertical direction and the area is 11.94 m2,as shown in Fig.19a.Furthermore,the entire grouting reinforced zone is a complete elliptical ring,with a long axis of 11.1 m in the vertical direction and a short axis of 10.8 m in the horizontal direction.However,when the lateral geostress increase to 4 MPa,the shape of the migration range of a single grouting hole changes to a circle with a radius of 1.4 m,as shown in Fig.19b.As σhfurther increases,as shown in Fig.19c,the shape of the migration range of a single grouting hole changes back to the ellipse again,but with a long axis of 1.75 m in the horizontal direction.Furthermore,some gaps appear between the isolated slurry migration islands at the left and right parts of the tunnels,while there are no gap appears at the top and bottom of the tunnel,where two continuous slurry migration zones are formed,respectively.Subsequently,as σhcontinuously increases from 8 to 12 MPa,the area of the migration zone in a single grouting hole continues to decrease from 4.12 to 2.90 m2,and a discrete slurry migration zone is formed,as shown in Fig.19d-f.

Fig.19.Slurry migration range in fractured strata models with different lateral ratios of geostress(σv:σh):(a)4 MPa:2 MPa,(b)4 MPa:4 MPa,(c)4 MPa:6 MPa,(d)4 MPa:8 MPa,(e) 4 MPa: 10 MPa,and (f) 4 MPa: 12 MPa.

As shown in Fig.20,it is clear that varying ratios of lateral geostress result in different grouting reinforcement effectiveness.For Case 7,under a relatively low lateral geostress (2 MPa),the convergence rate and the final convergence displacement of the surrounding rock mass have been effectively limited by the grouting reinforcement.At step 700,000,the area of the excavation section increases to 25.46 m2from 15.95 m2,compared with the tunnel without grouting reinforcement.When the lateral geostress increase to 4 MPa,horizontal and vertical stresses become equal,resulting in the formation of a complete circle ring grouting reinforcement zone.However,as illustrated in Table 12,the area of the grouting reinforced zone decreases from 92.79 m2to 69.95 m2,and the inner wall of the tunnel is partly reinforced by the injected slurry.This leads to the spalling of the broken rock mass at the upper and left part of the tunnel’s inner wall during the excavation process (as shown in Fig.20b).In Case 9,as the lateral geostress increase to 6 MPa,compared to the tunnel without grouting reinforcement,the area of the excavation face increases to 19.27 m2from 14.00 m2.Additionally,since the thickness of the annular reinforced zone in the horizontal direction is larger than the thickness in the vertical direction,the deformation at the left and right inner walls of the tunnel is well restrained and the convergence deformation mainly occurs at the top and bottom parts of the inner wall.In Cases 10 and 11,when the lateral geostress increases to 8 MPa and 10 MPa,the unreinforced gaps parallel to the horizontal direction begin to occur at the left and right parts of the tunnel.This behavior weakens the limitation on the convergence displacement on the lateral sides.Moreover,shear fractures begin to extend through the grouting-reinforced zone at the top and bottom of the tunnel,resulting in large deformation.In Case 12,where the horizontal stress increases to 12 MPa,not only does the crack in the grouting reinforcement zone at the top and bottom become more intensive,but also the shear cracks start to initiate in the grouting reinforcement zone at the two lateral sides of the tunnel.This leads to a more severe convergence of the tunnel section.At step 700,000,the area of the deformed section of the tunnels with and without grouting reinforcement is 13.1 m2and 11.03 m2,respectively,indicating that the effectiveness of grouting reinforcement is seriously weakened.

Table 12Slurry migration ranges of numerical Cases 7 to 12.

Fig.20.Comparison of convergence in tunnels with and without grouting reinforcement under different ratios of lateral geostress:(a)Case 7,(b)Case 8,(c)Case 9,(d)Case 10,(e)case 11,and (f) Case 12.

5.3.Influence of geostress levels on the grouting reinforcement effectiveness in tunnels with different fracture densities

All the models presented in the previous sections are conducted with the same fracture density.However,the grouting reinforcement effectiveness may vary depending on the geostress condition and the fracture density of fault fracture zones.Therefore,it is necessary to investigate the influence of geostress levels on the grouting reinforcement effectiveness in tunnels with different fracture densities.In this section,a series of excavation and grouting processes in the models with different fracture densities are conducted by NMM-TGR.Additionally,to characterize the fracture density of the stratum model in this study,the fracture density,Rf,of tunnel models is defined as the following equation:

whereNfis the number of fractures in the model,andStis the area of the model.According to Eq.(20),the fracture density of the tunnel model in the original manuscript is calculated and illustrated in Table 13.By multiplying the fracture density by a coefficient,Cf,three additional models with different fracture densities are obtained and are also illustrated in Table 13.Four tunnel models with varying fracture densities are subjected to the same boundary conditions(2 MPa and 16 MPa)and parameters as Cases 3 and 6 in Section 5.1.Based on the simulation results,the influence of geostress levels on the process of slurry migration and excavation in tunnels with different fracture densities is discussed below.

Table 13Tunnel models with various fracture densities and slurry migration ranges under different geostress.

The simulation results for the process of slurry migration are illustrated in Table 13 and Fig.21.When the fault fracture zone is with high fracture densities(Cf=100%),as shown in Fig.21a1 and b1,the connectivity of the slurry migration networks is high,and the slurry migration range is dominant by the geostress conditions.Therefore,under this circumstance,the area of the slurry migration zone decreases dramatically from 125.84 m2to 16.38 m2with an increase in geostress levels from 2 MPa to 16 MPa.Conversely,when the fault fracture zone is with low fracture densities(Cf=25%),as shown in Fig.21a4 and b4,the slurry migration range is mainly limited by the low connectivity of the slurry migration networks.Thus,under this circumstance,the area of the slurry migration zone remains consistent at 2.16 m2with an increase in geostress levels from 2 MPa to 16 MPa.

Fig.21.Slurry migration range in tunnel models with different fracture densities under different levels of geostress: (a) 2 MPa and (b) 16 MPa.

The simulation results for the excavation process of tunnels with and without grouting reinforcement are illustrated in Fig.22.For the models with a relatively high fracture density (Cf=100%),the slurry migration range decreases dramatically with an increase in geostress conditions.Therefore,the gaps between the two curves of the tunnel with and without grouting reinforcement,which indicates the effectiveness of grouting reinforcement,decrease dramatically with an increase in the geostress level (as shown in Fig.22a1 and b1).However,for the tunnel models with relatively low fracture densities(Cf=50% andCf=25%),the area of the slurry migration range is limited.As a result,the gap between the two curves is very small even under a low geostress condition,and the effectiveness of grouting reinforcement can be negligible under a high geostress condition.However,the variation in fracture density shows a significant influence on the deformation process of the tunnel excavation sections.As shown in Fig.22a and b,the area of the excavation section at step 700,000 increases with a decrease in the fracture density.This indicates that a decrease in fracture density can improve the integrity and bearing capacity of the surrounding rock mass.Thus,as shown in Fig.22,with a decrease in fracture density,the deformation rate of the surrounding rock mass in the tunnel with high geostress levels is slowed down,and the surrounding rock mass in the tunnel with low geostress levels can even remain undeformed.

Fig.22.Comparison of convergence in tunnel models with different fracture density under different levels of geostress: (a) 2 MPa (b) 16 MPa.

5.4.Discussion

Based on the numerical results presented above,it is evident that the effect of grouting reinforcement on broken rock masses in fault fracture zones with different fracture densities and geostress conditions varies significantly.Consequently,the optimization of the grouting reinforcement scheme should be tailored to the geostress conditions and fracture densities.Specifically,for the fault fracture zones with a relatively low level of geostress and a uniform fracture distribution,the slurry-reinforced zone formed is sufficient to cover the entire tunnel wall and prevent the large deformation of the surrounding rock mass.Therefore,in such cases,the optimization of the grouting scheme should focus on reducingAoand improvingReby reducing the number of grouting holes and increasing the spacing between them(as shown in Fig.23a and b).On the other hand,in the fault fracture zones with a high level of geostress and uniform fracture distribution,Amrsis relatively small and a grouting reinforced zone with many gaps is commonly encountered,resulting in a sharp decrease in grouting reinforcement effectiveness.Thus,in such cases,the optimization of the grouting scheme should focus on increasingAmrand eliminating the gaps between reinforced zones by increasing the number of grouting holes,reducing the spacing of grouting holes(as shown in Fig.23c and d) and increasing the slurry injection pressure.

Fig.23.Optimization of grouting hole arrangement: (a) reducing the number of grouting holes,(b) expanding the spacing of grouting holes,(c) adding the number of grouting holes,(d)reducing the spacing of grouting holes,(e)elliptical distribution and extra grouting holes at lateral sides,and(f)elliptical distribution with additional grouting holes at top and bottom.

For the fault fracture zones with different ratios of lateral geostress and a uniform fracture distribution,additional optimizations of the grouting scheme are required compared to the fault fracture zones under a hydrostatic pressure state.First,since the short axis of the elliptical grouting reinforced island is in the vertical direction,the unreinforced gaps parallel to the horizontal direction tend to occur at the lateral part of the tunnel (as shown in Fig.19c-f).Thus,to increase the continuity of the grouting reinforced zone,in addition to increasing the slurry pressures,it is also helpful to add more grouting holes at the lateral part of the tunnel in the circumferential direction as shown in Fig.23e.Moreover,as shown in Fig.19c-f,the shapes of the whole migration zone are no longer a circle but an ellipse with a short axis in the vertical direction,which indicates weaker reinforcement effectiveness in the tunnel top and bottom.Additionally,due to the main stress being in the horizontal direction,the crack initiation in the grouting reinforced zone no longer occurs uniformly around the tunnel and is more likely to occur at the top and bottom,resulting in larger convergence displacement at these locations (as shown in Fig.20c-f).Therefore,it is necessary to add more grouting holes at the top and bottom of the tunnel in the radial direction,as shown in Fig.23f,to better control crack initiation and tunnel deformation along the vertical direction.

For fault fracture zones with different fracture densities,it is necessary to implement different optimizations of the grouting scheme based on both the levels of geostress and fracture densities.Specifically,for the fault fracture zones with high fracture densities,the surrounding rock mass has low integrity and bearing capacity,and the slurry migration range is sensitive to the level of geostress.Therefore,when the high geostress condition is encountered,the key to optimizing the grouting scheme is to increaseAmrand eliminate the gaps between reinforced zones by adding more grouting holes,reducing the spacing of grouting holes(as shown in Fig.23c and d) and increasing the slurry injection pressures.Conversely,when the low geostress condition is encountered,since the area of slurry migration range at a single grouting hole is relatively large,the key to optimizing the grouting scheme is to reduceAoand improveReby reducing the number of grouting holes and increasing the spacing between grouting holes (as shown in Fig.23a and b).As for the fault fracture zone with low fracture densities,the surrounding rock mass has relatively high integrity and bearing capacity,and the slurry migration range is mainly limited by the low connectivity of fracture networks.Therefore,the advanced pre-grouting reinforcement adopted in this study is not recommended and the grouting reinforcement can be conducted in the surrounding rock mass behind the tunnel face,where the fracture networks are more developed and connected.

6.Conclusions

In this study,to investigate the influence of geostress conditions on the effectiveness of grouting reinforcement,an NMM-TGR method has been developed.The validity of the newly developed NMM-TGR method is established through a series of cases against analytical solutions as well as laboratory observations.Following the validation,the influence of geostress conditions on the effectiveness of grouting reinforcement is investigated by conducting a series of simulations about grouting reinforcement and tunnel excavation among the broken rock mass in fault fracture zones with various fracture densities under different geostress conditions.Based on the simulations results above,the following conclusions can be drawn:

(1) Due to the neglect of the influence of geostress condition,the slurry migration range is overestimated by 91% under a high level(16 MPa)of geostress condition.In addition,neglecting the influence of slurry pressure would also lead to an overestimation of the slurry migration range,and this overestimation can even reach 40% under low injecting pressures.

(2) For tunnels in the fault fracture zone with a quite low level of geostress and uniform fracture distribution,optimization of the grouting reinforcement scheme should focus on reducing the area of repeated grouting reinforced zone and cost of grouting by increasing the spacing between grouting holes.

(3) For tunnels in the fault fracture zone with a high level of geostress and uniform fracture distribution,the optimization of the grouting scheme should focus on increasing the area of the grouting reinforced zone and ensuring the continuity of the grouting reinforced zone by increasing the number of grouting holes,reducing the spacing of grouting holes,and increasing the slurry injecting pressures.

(4) For tunnels in the fault fracture zone with uniform fracture distribution and a high ratio of lateral geostress,extra optimization of the tunnel grouting scheme should be taken into consideration.First,extra grouting holes should be arranged on the two lateral sides of tunnels in the circumferential direction to eliminate the gaps between reinforced islands of each grouting hole.On the other hand,extra grouting holes should be arranged on the top and bottom of the tunnel in the radial direction to limit the crack initiation and deformation at the top and bottom of the tunnel.

(5) For tunnels in the fault fracture zone with high fracture density,the advance pre-grouting reinforcement in the surrounding rock mass ahead of the tunnel face is recommended.However,as for the tunnels in the fault fracture zone of low fracture density,it is recommended to conduct grouting reinforcement in the surrounding rock mass behind the tunnel face after the tunnel excavation.

Generally,this NMM-TGR method has demonstrated impressive potential in investigating the mechanism of the grouting reinforcement and the optimization of the grouting scheme.However,this paper is mainly focused on the influence of geostress conditions on the slurry migration,while its influence on the coagulation repair process,which is also influenced by geostress conditions,is not considered in the simulation process.In addition,since the deformation of broken surrounding rock mass cannot be well limited only by the approach of grouting reinforcement under a high level of geostress condition,other support approaches,such as lining supports and anchor supports,are usually required to be combined with the grouting approach.However,due to the limitation of 2D simulation,the compound supporting is not able to be simulated by this NMM-TGR method yet.Therefore,in the following studies,this NMM-TGR method will be further extended to 3D simulation for better representing the grouting reinforcement process and investigating the compound supporting mechanism between the grouting and other tunnel-supporting approaches.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the Guangdong Basic and Applied Basic Research Foundation (Grant No.2021A1515110304),the National Natural Science Foundation of China (Grant Nos.42077246 and 52278412).The authors are grateful for these financial supports.