Discontinuity development patterns and the challenges for 3D discrete fracture network modeling on complicated exposed rock surfaces

Natural slopes usually display complicated exposed rock surfaces that are characterized by complex and substantial terrain undulation and ubiquitous undesirable phenomena such as vegetation cover and rockfalls.This study presents a systematic outcrop research of fracture pattern variations in a complicated rock slope,and the qualitative and quantitative study of the complex phenomena impact on threedimensional(3D)discrete fracture network(DFN)modeling.As the studies of the outcrop fracture pattern have been so far focused on local variations,thus,we put forward a statistical analysis of global variations.The entire outcrop is partitioned into several subzones,and the subzone-scale variability of fracture geometric properties is analyzed(including the orientation,the density,and the trace length).The results reveal significant variations in fracture characteristics(such as the concentrative degree,the average orientation,the density,and the trace length)among different subzones.Moreover,the density of fracture sets,which is approximately parallel to the slope surface,exhibits a notably higher value compared to other fracture sets across all subzones.To improve the accuracy of the DFN modeling,the effects of three common phenomena resulting from vegetation and rockfalls are qualitatively analyzed and the corresponding quantitative data processing solutions are proposed.Subsequently,the 3D fracture geometric parameters are determined for different areas of the high-steep rock slope in terms of the subzone dimensions.The results show significant variations in the same set of 3D fracture parameters across different regions with density differing by up to tenfold and mean trace length exhibiting differences of 3e4 times.The study results present precise geological structural information,improve modeling accuracy,and provide practical solutions for addressing complex outcrop issues.

Algorithmic approach to discrete fracture network flow modeling in consideration of realistic connections in large-scale fracture networks

Analyzing rock mass seepage using the discrete fracture network(DFN)flow model poses challenges when dealing with complex fracture networks.This paper presents a novel DFN flow model that incorporates the actual connections of large-scale fractures.Notably,this model efficiently manages over 20,000 fractures without necessitating adjustments to the DFN geometry.All geometric analyses,such as identifying connected fractures,dividing the two-dimensional domain into closed loops,triangulating arbitrary loops,and refining triangular elements,are fully automated.The analysis processes are comprehensively introduced,and core algorithms,along with their pseudo-codes,are outlined and explained to assist readers in their programming endeavors.The accuracy of geometric analyses is validated through topological graphs representing the connection relationships between fractures.In practical application,the proposed model is employed to assess the water-sealing effectiveness of an underground storage cavern project.The analysis results indicate that the existing design scheme can effectively prevent the stored oil from leaking in the presence of both dense and sparse fractures.Furthermore,following extensive modification and optimization,the scale and precision of model computation suggest that the proposed model and developed codes can meet the requirements of engineering applications.

Development of an improved three-dimensional rough discrete fracture network model:Method and application

Structure plane is one of the important factors affecting the stability and failure mode of rock mass engineering.Rock mass structure characterization is the basic work of rock mechanics research and the important content of numerical simulation.A new 3-dimensional rough discrete fracture network(RDFN3D)model and its modeling method based on the Weierstrass-Mandelbrot(W-M)function were presented in this paper.The RDFN3D model,which improves and unifies the modelling methods for the complex structural planes,has been realized.The influence of fractal dimension,amplitude,and surface precision on the modeling parameters of RDFN3D was discussed.The reasonable W-M parameters suitable for the roughness coefficient of JRC were proposed,and the relationship between the mathematical model and the joint characterization was established.The RDFN3D together with the smooth 3-dimensional discrete fracture network(DFN3D)models were successfully exported to the drawing exchange format,which will provide a wide application in numerous numerical simulation codes including both the continuous and discontinuous methods.The numerical models were discussed using the COMSOL Multiphysics code and the 3-dimensional particle flow code,respectively.The reliability of the RDFN3D model was preliminarily discussed and analyzed.The roughness and spatial connectivity of the fracture networks have a dominant effect on the fluid flow patterns.The research results can provide a new geological model and analysis model for numerical simulation and engineering analysis of jointed rock mass.

Borehole stability in naturally fractured rocks with drilling mud intrusion and associated fracture strength weakening:A coupled DFN-DEM approach

Borehole instability in naturally fractured rocks poses significant challenges to drilling.Drilling mud invades the surrounding formations through natural fractures under the difference between the wellbore pressure(P w)and pore pressure(P p)during drilling,which may cause wellbore instability.However,the weakening of fracture strength due to mud intrusion is not considered in most existing borehole stability analyses,which may yield significant errors and misleading predictions.In addition,only limited factors were analyzed,and the fracture distribution was oversimplified.In this paper,the impacts of mud intrusion and associated fracture strength weakening on borehole stability in fractured rocks under both isotropic and anisotropic stress states are investigated using a coupled DEM(distinct element method)and DFN(discrete fracture network)method.It provides estimates of the effect of fracture strength weakening,wellbore pressure,in situ stresses,and sealing efficiency on borehole stability.The results show that mud intrusion and weakening of fracture strength can damage the borehole.This is demonstrated by the large displacement around the borehole,shear displacement on natural fractures,and the generation of fracture at shear limit.Mud intrusion reduces the shear strength of the fracture surface and leads to shear failure,which explains that the increase in mud weight may worsen borehole stability during overbalanced drilling in fractured formations.A higher in situ stress anisotropy exerts a significant influence on the mechanism of shear failure distribution around the wellbore.Moreover,the effect of sealing natural fractures on maintaining borehole stability is verified in this study,and the increase in sealing efficiency reduces the radial invasion distance of drilling mud.This study provides a directly quantitative prediction method of borehole instability in naturally fractured formations,which can consider the discrete fracture network,mud intrusion,and associated weakening of fracture strength.The information provided by the numerical approach(e.g.displacement around the borehole,shear displacement on fracture,and fracture at shear limit)is helpful for managing wellbore stability and designing wellbore-strengthening operations.

Characterizing the influence of stress-induced microcracks on the laboratory strength and fracture development in brittle rocks using a finite-discrete element method-micro discrete fracture network FDEM-μDFN approach

Heterogeneity is an inherent component of rock and may be present in different forms including mineralheterogeneity, geometrical heterogeneity, weak grain boundaries and micro-defects. Microcracks areusually observed in crystalline rocks in two forms： natural and stress-induced; the amount of stressinducedmicrocracking increases with depth and in-situ stress. Laboratory results indicate that thephysical properties of rocks such as strength, deformability, P-wave velocity and permeability areinfluenced by increase in microcrack intensity. In this study, the finite-discrete element method （FDEM）is used to model microcrack heterogeneity by introducing into a model sample sets of microcracks usingthe proposed micro discrete fracture network （mDFN） approach. The characteristics of the microcracksrequired to create mDFN models are obtained through image analyses of thin sections of Lac du Bonnetgranite adopted from published literature. A suite of two-dimensional laboratory tests including uniaxial,triaxial compression and Brazilian tests is simulated and the results are compared with laboratory data.The FDEM-mDFN models indicate that micro-heterogeneity has a profound influence on both the mechanicalbehavior and resultant fracture pattern. An increase in the microcrack intensity leads to areduction in the strength of the sample and changes the character of the rock strength envelope. Spallingand axial splitting dominate the failure mode at low confinement while shear failure is the dominantfailure mode at high confinement. Numerical results from simulated compression tests show thatmicrocracking reduces the cohesive component of strength alone, and the frictional strength componentremains unaffected. Results from simulated Brazilian tests show that the tensile strength is influenced bythe presence of microcracks, with a reduction in tensile strength as microcrack intensity increases. Theimportance of microcrack heterogeneity in reproducing a bi-linear or S-shape failure envelope and itseffects on the mechanisms leading to spalling damage near an underground opening are also discussed.

Fractured reservoir modeling by discrete fracture network and seismic modeling in the Tarim Basin,China

Fractured reservoirs are an important target for oil and gas exploration in the Tarim Basin and the prediction of this type of reservoir is challenging.Due to the complicated fracture system in the Tarim Basin,the conventional AVO inversion method based on HTI theory to predict fracture development will result in some errors.Thus,an integrated research concept for fractured reservoir prediction is put forward in this paper.Seismic modeling plays a bridging role in this concept,and the establishment of an anisotropic fracture model by Discrete Fracture Network （DFN） is the key part.Because the fracture system in the Tarim Basin shows complex anisotropic characteristics,it is vital to build an effective anisotropic model.Based on geological,well logging and seismic data,an effective anisotropic model of complex fracture systems can be set up with the DFN method.The effective elastic coefficients,and the input data for seismic modeling can be calculated.Then seismic modeling based on this model is performed,and the seismic response characteristics are analyzed.The modeling results can be used in the following AVO inversion for fracture detection.

Application of laser scanning for rock mass characterization and discrete fracture network generation in an underground limestone mine

Terrestrial laser scanning(TLS) is a useful technology for rock mass characterization. A laser scanner produces a massive point cloud of a scanned area, such as an exposed rock surface in an underground tunnel,with millimeter precision. The density of the point cloud depends on several parameters from both the TLS operational conditions and the specifications of the project, such as the resolution and the quality of the laser scan, the section of the tunnel, the distance between scanning stations, and the purpose of the scans. One purpose of the scan can be to characterize the rock mass and statistically analyze the discontinuities that compose it for further discontinuous modeling. In these instances, additional data processing and a detailed analysis should be performed on the point cloud to extract the parameters to define a discrete fracture network(DFN) for each discontinuity set. I-site studio is a point cloud processing software that allows users to edit and process laser scans. This software contains a set of geotechnical analysis tools that assist engineers during the structural mapping process, allowing for greater and more representative data regarding the structural information of the rock mass, which may be used for generating DFNs. This paper presents the procedures used during a laser scan for characterizing discontinuities in an underground limestone mine and the results of the scan as applied to the generation of DFNs for further discontinuous modeling.

Identifying Non-Darcian Flow and Non-Fickian Pressure Propagation in Field-Scale Discrete Fracture Networks

Non-Darcian flow has been well documented for fractured media, while the potential non-Darcian flow and its driven factors in field-scale discrete fracture networks (DFNs) remain obscure. This study conducts Monte Carlo simulations of water flow through DFNs to identify non-Darcian flow and non-Fickian pressure propagation in field-scale DFNs, by adjusting fracture density, matrix hydraulic conductivity, and the general hydraulic gradient. Numerical simulations and analyses show that interactions of the fracture architecture with the hydraulic gradient affect non-Darcian flow in DFNs, by generating and adjusting complex pathways for water. The fracture density affects significantly the propagation of hydraulic head/pressure in the DFN, likely due to fracture connectivity and flow channeling. The non-Darcian flow pattern may not be directly correlated to the non-Fickian pressure propagation process in the regional-scale DFNs, because they refer to different states of water flow and their controlling factors may not be the same. Findings of this study improve our understanding of the nature of flow in DFNs.

Well interference evaluation considering complex fracture networks through pressure and rate transient analysis in unconventional reservoirs

Severe well interference through complex fracture networks(CFNs)can be observed among multi-well pads in low permeability reservoirs.The well interference analysis between multi-fractured horizontal wells(MFHWs)is vitally important for reservoir effective development.Well interference has been historically investigated by pressure transient analysis,while it has shown that rate transient analysis has great potential in well interference diagnosis.However,the impact of complex fracture networks(CFNs)on rate transient behavior of parent well and child well in unconventional reservoirs is still not clear.To further investigate,this paper develops an integrated approach combining pressure and rate transient analysis for well interference diagnosis considering CFNs.To perform multi-well simulation considering CFNs,non-intrusive embedded discrete fracture model approach was applied for coupling fracture with reservoir models.The impact of CFN including natural fractures and frac-hits on pressure and rate transient behavior in multi-well system was investigated.On a logelog plot,interference flow and compound linear flow are two new flow regimes caused by nearby producers.When both NFs and frac-hits are present in the reservoir,frac-hits have a greater impact on well#1 which contains frac-hits,and NFs have greater impact on well#3 which does not have frac-hits.For all well producing circumstances,it might be challenging to see divergence during pseudosteady state flow brought on by frac-hits on the logelog plot.Besides,when NFs occur,reservoir depletion becomes noticeable in comparison to frac-hits in pressure distribution.Application of this integrated approach demonstrates that it works well to characterize the well interference among different multi-fractured horizontal wells in a well pad.Better reservoir evaluation can be acquired based on the new features observed in the novel model,demonstrating the practicability of the proposed approach.The findings of this study can help for better evaluating well interference degree in multi-well systems combing PTA and RTA,which can reduce the uncertainty and improve the accuracy of the well interference analysis based on both field pressure and rate data.

Horizontal Well Interference Performance and Water Injection Huff and Puff Effect on Well Groups with Complex Fracture Networks:A Numerical Study

Well interference has become a common phenomenon with the increasing scale of horizontal well fracturing.Recent studies on well interference in horizontal wells do not properly reflect the physical model of the postfracturing well groups and the realistic fracturing process of infill wells.Establishing the correspondence between well interference causative factors and manifestations is of great significance for infill well deployment and secondary oil recovery.In this work,we develop a numerical model that considers low velocity non-Darcy seepage inshale reservoirs to study the inter-well interferencephenomenon that occurs in theSantanghufield,andconstruct an explicit hydraulic fracture and complex natural fracture network model with an embedded discrete fracture model,focusing on the effect of fracture network morphology on well interactions.The model also considers a multi-segment wellbore model to accommodate the effect of inter-well crossflow on wellbore tubular flow.The changes in formation pressure and water saturation during fracturing are performed by controlling the injection pressure and water injection rate.The result shows that the shape of the fracture network generated by the infill well with the old well determines the subsequent fluid and oil-increasing performance of the disturbed well.The synergistic production or competitive relationship formed by fractures with different connectivity between the two wells determines the positive and negative effects of the interference.The paper also investigates the adaptation study of water injection huff and puff schemes for well groups with different connectivity,and demonstrated a potential yield increase of up to 10.85%under adaptation injection.This method of identifying well interference based on the production dynamics of affected wells and the subsequent corresponding water injection method provides valuable references for the selection of secondary oil recovery measures.

A discrete model for prediction of radon flux from fractured rocks

Prediction of radon flux from the fractured zone of a propagating cave mine is basically associated with uncertainty and complexity. For instance, there is restricted access to these zones for field measure- ments, and it is quite difficult to replicate the complex nature of both natural and induced fractures in these zones in laboratory studies. Hence, a technique for predicting radon flux from a fractured rock using a discrete fracture network （DFN） model is developed to address these difficulties. This model quantifies the contribution of fractures to the total radon flux, and estimates the fracture density from a measured radon flux considering the effects of advection, diffusion, as well as radon generation and decay. Radon generation and decay are classified as reaction processes. Therefore, the equation solved is termed as the advection-diffusion-reaction equation （ADRE）. Peclet number （Pe）, a conventional dimensionless parameter that indicates the ratio of mass transport by advection to diffusion, is used to classify the transport regimes. The results show that the proposed model effectively predicts radon flux from a fractured rock. An increase in fracture density for a rock sample with uniformly distributed radon generation rate can elevate radon flux significantly compared with another rock sample with an equivalent increase in radon generation rate. In addition to Pe, two other independent dimensionless parameters （derived for radon transport through fractures） significantly affect radon dimensionless flux. Findings provide insight into radon transport through fractured rocks and can be used to improve radon control measures for proactive mitigation.

Evaluating the mechanical properties of anisotropic shale containing bedding and natural fractures with discrete element modeling

Natural fracture data from one of the Carboniferous shale masses in the eastern Qaidam Basin were used to establish a stochastic model of a discrete fracture network and to perform discrete element simulation research on the size efect and mechanical parameters of shale.Analytical solutions of fctitious joints in transversely isotropic media were derived,which made it possible for the proposed numerical model to simulate the bedding and natural fractures in shale masses.The results indicate that there are two main factors infuencing the representative elementary volume(REV)size of a shale mass.The frst and most decisive factor is the presence of natural fractures in the block itself.The second is the anisotropy ratio:the greater the anisotropy is,the larger the REV.The bedding angle has little infuence on the REV size,whereas it has a certain infuence on the mechanical parameters of the rock mass.When the bedding angle approaches the average orientation of the natural fractures,the mechanical parameters of the shale blocks decrease greatly.The REV representing the mechanical properties of the Carboniferous shale masses in the eastern Qaidam Basin were comprehensively identifed by considering the infuence of bedding and natural fractures.When the numerical model size is larger than the REV,the fractured rock mass discontinuities can be transformed into equivalent continuities,which provides a method for simulating shale with natural fractures and bedding to analyze the stability of a borehole wall in shale.

Numerical investigation of dual-porosity model with transient transfer function based on discrete-fracture model

Based on the characteristics of fractures in naturally fractured reservoir and a discrete-fracture model, a fracture network numerical well test model is developed. Bottom hole pressure response curves and the pressure field are obtained by solving the model equations with the finite-element method. By analyzing bottom hole pressure curves and the fluid flow in the pressure field, seven flow stages can be recognized on the curves. An upscaling method is developed to compare with the dual-porosity model （DPM）. The comparisons results show that the DPM overestimates the inter-porosity coefficient ）, and the storage factor w. The analysis results show that fracture conductivity plays a leading role in the fluid flow. Matrix permeability influences the beginning time of flow from the matrix to fractures. Fractures density is another important parameter controlling the flow. The fracture linear flow is hidden under the large fracture density. The pressure propagation is slower in the direction of larger fracture density.

Estimation of fracture size and azimuth in the universal elliptical disc model based on trace information

The geometric characteristics of fractures within a rock mass can be inferred by the data sampling from boreholes or exposed surfaces.Recently,the universal elliptical disc(UED)model was developed to represent natural fractures,where the fracture is assumed to be an elliptical disc and the fracture orientation,rotation angle,length of the long axis and ratio of short-long axis lengths are considered as variables.This paper aims to estimate the fracture size-and azimuth-related parameters in the UED model based on the trace information from sampling windows.The stereological relationship between the trace length,size-and azimuth-related parameters of the UED model was established,and the formulae of the mean value and standard deviation of trace length were proposed.The proposed formulae were validated via the Monte Carlo simulations with less than 5%of error rate between the calculated and true values.With respect to the estimation of the size-and azimuth-related parameters using the trace length,an optimization method was developed based on the pre-assumed size and azimuth distribution forms.A hypothetical case study was designed to illustrate and verify the parameter estimation method,where three combinations of the sampling windows were used to estimate the parameters,and the results showed that the estimated values could agree well with the true values.Furthermore,a hypothetical three-dimensional(3D)elliptical fracture network was constructed,and the circular disc,non-UED and UED models were used to represent it.The simulated trace information from different models was compared,and the results clearly illustrated the superiority of the proposed UED model over the existing circular disc and non-UED models。

Analysis of fracture propagation and shale gas production by intensive volume fracturing

This paper presents an integrated study from fracture propagation modeling to gas flow modeling and a correlation analysis to explore the key controlling factors of intensive volume fracturing.The fracture propagation model takes into account the interaction between hydraulic fracture and natural fracture by means of the displacement discontinuity method(DDM)and the Picard iterative method.The shale gas flow considers multiple transport mechanisms,and the flow in the fracture network is handled by the embedded discrete fracture model(EDFM).A series of numerical simulations are conducted to analyze the effects of the cluster number,stage spacing,stress difference coefficient,and natural fracture distribution on the stimulated fracture area,fractal dimension,and cumulative gas production,and their correlation coefficients are obtained.The results show that the most influential factors to the stimulated fracture area are the stress difference ratio,stage spacing,and natural fracture density,while those to the cumulative gas production are the stress difference ratio,natural fracture density,and cluster number.This indicates that the stress condition dominates the gas production,and employing intensive volume fracturing(by properly increasing the cluster number)is beneficial for improving the final cumulative gas production.

Assessing fracturing mechanisms and evolution of excavation damaged zone of tunnels in interlocked rock masses at high stresses using a finitediscrete element approach

Deep underground excavations within hard rocks can result in damage to the surrounding rock mass mostly due to redistribution of stresses.Especially within rock masses with non-persistent joints,the role of the pre-existing joints in the damage evolution around the underground opening is of critical importance as they govern the fracturing mechanisms and influence the brittle responses of these hard rock masses under highly anisotropic in situ stresses.In this study,the main focus is the impact of joint network geometry,joint strength and applied field stresses on the rock mass behaviours and the evolution of excavation induced damage due to the loss of confinement as a tunnel face advances.Analysis of such a phenomenon was conducted using the finite-discrete element method (FDEM).The numerical model is initially calibrated in order to match the behaviour of the fracture-free,massive Lac du Bonnet granite during the excavation of the Underground Research Laboratory (URL) Test Tunnel,Canada.The influence of the pre-existing joints on the rock mass response during excavation is investigated by integrating discrete fracture networks (DFNs) of various characteristics into the numerical models under varying in situ stresses.The numerical results obtained highlight the significance of the pre-existing joints on the reduction of in situ rock mass strength and its capacity for extension with both factors controlling the brittle response of the material.Furthermore,the impact of spatial distribution of natural joints on the stability of an underground excavation is discussed,as well as the potentially minor influence of joint strength on the stress induced damage within joint systems of a non-persistent nature under specific conditions.Additionally,the in situ stress-joint network interaction is examined,revealing the complex fracturing mechanisms that may lead to uncontrolled fracture propagation that compromises the overall stability of an underground excavation.

Numerical Modeling of Natural Fracture Distributions in Shale

The production efficiency of shale gas is affected by the interaction between hydraulic and natural fractures.This study presents a simulation of natural fractures in shale reservoirs,based on a discrete fracture network(DFN)method for hydraulic fracturing engineering.Fracture properties of the model are calculated from core fracture data,according to statistical mathematical analysis.The calculation results make full use of the quantitative information of core fracture orientation,density,opening and length,which constitute the direct and extensive data of mining engineering.The reliability and applicability of the model are analyzed with regard to model size and density,a calculation method for dominant size and density being proposed.Then,finite element analysis is applied to a hydraulic fracturing numerical simulation of a shale fractured reservoir in southeastern Chongqing.The hydraulic pressure distribution,fracture propagation,acoustic emission information and in situ stress changes during fracturing are analyzed.The results show the application of fracture statistics in fracture modeling and the influence of fracture distribution on hydraulic fracturing engineering.The present analysis may provide a reference for shale gas exploitation.

Multi-dimensional size effects and representative elements for nonpersistent fractured rock masses: A perspective of geometric parameter distribution

This study takes a fractured rock mass in the Datengxia Hydropower Station,China as an example to analyze the size effects and determine the representative elementary sizes.A novel method considering geometric parameter distributions is proposed in this work.The proposed method can quickly and simply determine the size effects and representative elementary sizes.Specifically,geometric parameter distributions,including fracture frequency,size and orientation,are generated on the basis of the Bernoulli trial and Monte Carlo simulation.The distributions are assessed using the coefficient of variation(CV),and the acceptable variations for CV(5%,10%and 20%)are used to determine representative elementary sizes.Generally,the representative element of rock masses is the representative elementary volume(REV).The present study extends the representative element to other dimensions,i.e.representative elementary length(REL)and representative elementary area(REA)for one and two dimensions,respectively.REL and REA are useful in studying the size effects of one-(1D)and twodimensional(2D)characteristics of rock masses.The relationships among multi-dimensional representative elementary sizes are established.The representative elementary sizes reduce with the increase in the dimensions,and REA and REV can be deduced by REL.Therefore,the proposed method can quickly and simply determine REL and further estimate REA and REV,which considerably improves the efficiency of rock mass analysis.

Multi-domain equivalent method for prediction of elastic modulus of complex fractured rock mass

Reliable estimation of deformation and failure behaviors of fractured rock mass is important for practical engineering design.This study proposes a multi-domain equivalent method for fracture network to estimate the deformation properties of complex fractured rock mass.It comprehends both the advantages of the discrete fracture network model and the equivalent continuum model to capture the features of discontinuities explicitly while reducing computational intensity.The complex fracture network is stochastically split into a number of subfracture networks according to the domain,length or angle.An analytical solution is derived to infer theoretically the relationship between the elastic moduli of the original complex fractured rock mass and the split subfractured rock masses by introducing a correction term based on the deformation superposition principle.Numerical simulations are conducted to determine the elastic moduli of split subfractured rock masses using universal distinct element code(UDEC),while the elastic modulus of the original model is estimated based on the currently proposed analytical relationship.The results show that the estimation accuracy with the current domainbased splitting model is far superior compared to those with the other two splitting models.Thus,the estimation method of elastic modulus of complex fractured rock mass based on domain splitting mode of fracture network is identified as the multi-domain equivalent method proposed in this paper.The reliability of this method is evaluated,and its high computational efficiency is demonstrated through exemplification with regard to different geometric configurations for stochastically artificial discrete fracture network.The proposed multi-domain equivalent method constructs the theoretical framework except for the regression analysis hypothesis compared to the density-reduced model equivalent method.

Complex fracture propagation model and plugging timing optimization for temporary plugging fracturing in naturally fractured shale

In this paper,a viscoelasticity-plastic damage constitutive equation for naturally fractured shale is deduced,coupling nonlinear tensile-shear mixed fracture mode.Dynamic perforation-erosion on fluid re-distribution among multi-clusters are considered as well.DFN-FEM(discrete fracture network combined with finite element method)was developed to simulate the multi-cluster complex fractures propagation within temporary plugging fracturing(TPF).Numerical results are matched with field injection and micro-seismic monitoring data.Based on geomechanical characteristics of Weiyuan deep shale gas reservoir in Sichuan Basin,SW China,a multi-cluster complex fractures propagation model is built for TPF.To study complex fractures propagation and the permeability-enhanced region evolution,intersecting and competition mechanisms between the fractures before and after TPF treatment are revealed.Simulation results show that:fracture from middle cluster is restricted by the fractures from side-clusters,and side-clusters plugging is benefit for multi fractures propagation in uniformity;optimized TPF timing should be delayed within a higher density or strike of natural fractures;Within a reservoir-featured natural fractures distribution,optimized TPF timing for most clustered method is 2/3 of total fluid injection time as the optimal plugging time under different clustering modes.