Xiao Chen

Chong Shi

Yun Jia

Yi-Ping Zhang

Huai-Ning Ruan

Numerical study on seepage properties of rock mass with non-penetrating fracture using discrete element method

Keywords: Discrete element method, Non-penetrating fracture, Pipe flow algorithm, Pore pressure, Permeability, Fractured rock mass

This paper aims to investigate the seepage process of fractured rock mass at a mesoscopic level, using the discrete element method. To this end, an improved pipe flow algorithm is proposed within the framework of particle flow code, including the heterogeneity of initial micro-cracks and the deformation of solid skeleton. A series of numerical samples with a single non-penetrating fracture are then studied. The relationships between the rock permeability and different properties of fracture (i.e. the aperture, dip angle and length) are investigated. The influence of fracture properties on the distribution of pore pressure in the rock mass is also analyzed. The obtained numerical results exhibit that the seepage properties of rock masse with non-penetrating fracture depend strongly on the dip angle and length of fracture and are slightly disturbed by the fracture aperture. Moreover, its permeability is almost of the same order as that of non-fractured rock. The numerical results and mesoscale analysis can help us understand the influence of fracture properties on the permeability of host rock and provide a practical guidance for the prevention of hydraulic disasters and the 2 development of mineral resources.

Introduction

In geotechnical engineering, the presence of various discontinuities significantly influences the seepage characteristics of rock mass. For instance, in petroleum engineering, the creation and propagation of fractures can enhance the efficiency of oil and gas exploration, but they may also pose risks to underground structures due to gushing groundwater. Therefore, understanding the influence of fracture properties on the seepage behavior of fractured rock mass is crucial for the exploration of natural resources and the safety of underground structures.

In the literature, a number of laboratory experiments have been carried out to study the seepage characteristics of rocks (Heiland 2003a;[START_REF] Wang | Numerical Modeling of Pore Pressure Influence on Fracture Evolution in Brittle Heterogeneous Rocks Rock Mechanics and Rock Engineering[END_REF]. Most of these studies focus on the permeability and failure characteristics of rocks during the loading and unloading process (Heiland 2003b). Recently, [START_REF] Xu | Permeability evolution of sandstone under short-term and long-term triaxial compression[END_REF] measured the short-term and long-term seepage characteristics of rock under different confining pressures using the transient method. They found that the evolution of permeability under different confining and osmotic pressures was similar, and a linear relationship could be established between the rock permeability and its volumetric strain. Liu et al. (2020a) analyzed the permeability evolution of some sandstones, including the dependence of permeability on the microstructure. However, due to the microstructural complexity of rock masses, it's difficult to quantify the influence of fracture on the seepage behavior of rock mass. In view of this, some scholars studied the permeability characteristics of a single fracture. [START_REF] Qian | Experimental evidence of scale-dependent hydraulic conductivity for fully developed turbulent flow in a single fracture[END_REF] investigated the hydraulic conductivity of a single fracture with various openings and roughness under different hydraulic gradients. Zhang et al. (2013a) investigated the effects of fluid viscosity, fracture opening, fracture inclination, and driving velocity on the relative permeability of a single fracture. Recently, [START_REF] Wu | Investigations on the seepage characteristics of fractured sandstone based on NMR real-time imaging[END_REF] carried out a series of seepage experiments based on some sandstone samples with different fracture characteristics, monitored the real-time fracture seepage process, and discussed the relationship between the inclination angle and number of fractures and permeability. However, the fractured samples in the above-mentioned laboratory tests are treated as two independent entities (i.e. fracture and rock matrix) and the fractures used in these studies are generally penetrating. In practice, most of fractures in rock mass are internal defects and embedded, rather than penetrating.

Considering the insufficiency of laboratory experiments, some numerical methods have been used to study the seepage characteristics of fractured rock masses. Most of these methods are based on the concept of continuum mechanics, such as the embedded finite element method-EFEM [START_REF] Huang | Experimental and Numerical Studies of the Hydraulic Properties of Three-Dimensional Fracture Networks with Spatially Distributed Apertures[END_REF], finite difference method-FDM [START_REF] Chen | Numerical study on the effects of fracture parameters on permeability in fractured rock with extremely low matrix permeability[END_REF]) and extended finite element method-XFEM [START_REF] Dehghan | 3D Numerical Modeling of the Propagation of Hydraulic Fracture at Its Intersection with Natural (Pre-existing)[END_REF], etc. However, these approaches commonly struggle to predict the fracture propagation and interaction. In parallel, some methods based on discontinuity mechanics have also been used, such as numerical manifold method-NMM [START_REF] Wu | Hydro-mechanical simulation of the saturated and semi-saturated porous soil-rock mixtures using the numerical manifold method Computer[END_REF][START_REF] Yang | Modeling the entire progressive failure process of rock slopes using a strength-based criterion Computers and[END_REF].

The NMM has been applied to simulate two-dimensional fluid driven fracturing in rocks by Yang et al. (2018) and analyze the evolution of three-dimensional fractures in concrete structures [START_REF] Yang | Three-dimensional fracture propagation with numerical manifold method Engineering Analysis with Boundary[END_REF]. Although the NMM combines the advantages of finite element method and discontinuity mechanics, its application in practical engineering is difficult due to the high complexity of its algorithm. On the other hand, discrete element method (DEM), due to its explicit numerical scheme and simple motion equations, has been largely used to study the mechanical properties and microscopic failure mechanism of geomaterials (Cundall and Hart 1993;[START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Hu | Parametric study of the smooth-joint contact model on the mechanical behavior of jointed rock International Journal for Numerical and[END_REF]. While different models have been proposed to study the fluid-structure interaction in the framework of discrete element method, most of them require the total integration of particle and fluid calculation, as well as the transition of particle information to fluid schema. To overcome these difficulties, the pipe network flow model has been implemented in the particle flow code (PFC). Al-Busaidi et al. (2005) have succussefully used the pipe network flow model to reproduce the fluid flow injected into a low-permeability porous media. It has also been used to study hydraulic fracturing problems, including the dependence of seepage process on the fluid viscosity, fluid injection rate, insitu stress state and distribution of particle size (Liu et al. 2020b;[START_REF] Shimizu | The distinct element analysis for hydraulic fracturing in hard rock considering fluid viscosity and particle size distribution[END_REF][START_REF] Wang | Simulation of hydraulic fracturing using particle flow method and application in a coal mine[END_REF]Zhang et al. 2013b;[START_REF] Zhou | Numerical investigation of fluid-driven near-borehole fracture propagation in laminated reservoir rock using PFC 2[END_REF]. The interaction modes of natural fractures and hydraulic fractures under different boundary conditions were also studied by using the pipe network flow model [START_REF] Wang | The effect of natural fractures on hydraulic fracturing propagation in coal seams[END_REF][START_REF] Zhao | Controls of natural fractures on the texture of hydraulic fractures in rock[END_REF]. However, the original fluid-solid coupling algorithm embedded in PFC is questionable, leading to the proposal of several improved algorithm. For instance, some modified algorithms considering the failure of contact bond have been proposed. To reproduce the pressure evolution related to the crack propagation, [START_REF] Zhang | An improved hydromechanical model for particle flow simulation of fractures in fluid-saturated rocks[END_REF] recently introduced a new update equation for hydraulic fracture aperture in the hydro-mechanical coupling process. However, most of these research works focus on the hydraulic fracturing issues. When the pipe flow model is applied in hydraulic fracturing problems, pore pressure undergoes continuous fluctuations, eliminating the need for a steady-state condition. In contrast, when simulating the seepage behavior of rock masses, the steady-state seepage conditions are required. This necessitates the utilization of more precise parameters and intricate algorithmic details, significantly elevating the complexity of the algorithm.

Although some numerical studies have been conducted on the permeability of fractured rock mass, most of them focus on the fracture network, and the effects of single embedded fracture have not been studied. Moreover, the fractures studied in the previous research are penetrating, and rare research have been performed on the non-penetrating embedded fracture on the seepage processes of rock mass. To the best knowledge of authors, the pipe flow model has been widely used in hydraulic fracturing, but has rarely been used to study the seepage behavior of rock mass. In the present work, the pipe flow model will be used to numerically investigate the seepage behavior of rock masses with single non-penetrating embedded fractures. Furthermore, an improved pipe flow algorithm will be proposed to take into account the pressure evolution related to the deformation of solid skeleton.

Additionally, the analysis of the effect of single non-penetrating embedded fractures on the seepage behavior of rock masses, performed in the present work, is capable of providing novel insights into the field of fractured rock mechanics.

In the present work, the seepage characteristics of rock mass with a single non-penetrating embedded fracture is numerically investigated, with special attention to the fracture properties. To accomplish this, an improved pipe flow algorithm is proposed and implemented in the discrete element modeling framework of PFC2D. The single-fracture rock samples with different fracture properties are then investigated. Based on the numerical results of permeability and pore pressure distributions, the effects of embedded fracture on the seepage properties of rock mass are discussed.

Modeling methods and theories

In the particle discrete element method, the intact rock is presented as a dense packing of nondeformable spheres, bonded with contacts. The bonding state of a geotechnical material is represented by the contacts between spheres. As a result, the macroscopic mechanical properties of different geotechnical materials can be obtained by applying different meso-parameters on the spheres and contacts.

Initially, contacts are established at the known positions of particles. Subsequently, the forcedisplacement law is utilized for each contact to update the contact forces between the particles [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Fan | Cone complimentary-based numerical manifold method modeling frictional and cohesive contact problems[END_REF][START_REF] Wu | Enriched three-field numerical manifold formulation for dynamics of fractured saturated porous media[END_REF]. Once the resultant force and moment of the particle are determined, the particle's velocity and displacement can be obtained using the motion law. The equation governing translational motion can be expressed as follows:

(1) where m is the particle mass, and represents the particle acceleration. F indicates the sum of the applied forces acting on the particle, including the body force Fb, the applied external force Fa and the contact force Fc. In practice, Fb and Fa are generally treated as constants while Fc needs to be recalculated at each time step based on the constitutive model of contact. The linear parallel bond model, commonly used in the simulation of mechanical behavior of rock-like materials, is adopted in the present work. As shown in Fig. 1, the force-displacement law for the parallel bond model can described as follows:

,

,

where Fn and Fs are the normal and tangential contact forces, respectively. kn denotes the normal stiffness and ks represents the shear stiffness. Δδn and Δδs represent the increments of normal and shear displacement at each timestep, respectively. A refers to the area of the contact plane between two particles. Then, the tensile stress(σ) and shear stress(τ) of the bond element are calculated as follows:

,

,

where I and Mb are the moment of inertia and the bending moment of the contact plain, respectively. R is the average radius of two bonded particles, and β is the moment contribution coefficient.

The bond breaks in tension when the contact stress exceeds the tensile strength:

σ > σc (6)
where σc is the tensile-strength limit of the bonded element. If the bond has not broken in tension, the shear-strength limit is enforced. Once the shear stress exceeds the shear-strength limit, the bond breaks in shear mode :

τc=c-σtanφ

(7) n = + n n n F F k A d D s s s s F F k A d = -D = + n b F M R A I s b s F A t =
where τc is the shear-strength limit of the bonded element, c and φ are the cohesion and friction angle of the contact element.

Review of pipe flow model

In geotechnical materials, such as rock, fluid flow occurs in two types of voids, called "pores" and "throats". Pores refer to larger empty spaces that are surrounded by multiple skeleton particles, making up the material skeleton. These pores primarily contribute to fluid storage within the material.

On the other hand, throats are narrow passageways located between two particles or the constricted spaces between adjacent pores. These throats play a crucial role in facilitating the movement of fluids through the material.

In the particle discrete model, a large number of pores exist between the particles in the material assembly (Fig. 2). These pores are considered as domains within the pipe network flow model, representing the storage capacity for fluids. Since the particles in the assembly are generally in close contact with each other, identifying individual throats becomes challenging in the numerical simulation. To overcome this challenge, a numerical approach employs pipes to represent throats between two particles. By connecting the particles with pipes, a network of fluid flow channels is formed, allowing for the simulation of fluid transport through the material (Fig. 3a). This simplification is capable of modeling the fluid flow behavior with the consideration of the complex geometry of the material. Under two-dimensional conditions, the fluid rate q in the pipe can be calculated as follows:

(8)

where ɑ and Lp represent the aperture and length of the pipe, respectively. μ refers to the fluid viscosity, and ΔP denotes the pressure difference between two adjacent domains. The pipe length Lp is determined by the sum of the radii of two adjacent particles on both sides of the pipe. Regarding the pipe aperture ɑ, when the sample is free to external loading, it takes the initial value, denoted as ɑ0.

After that, as the applied normal force Fn is compressible, the fracture aperture decreases asymptotically toward zero with the increase of Fn. This means that as the normal force becomes larger, the size of the fracture aperture diminishes, resulting in a decrease in the fluid flow capacity Where represents the compressive normal force when the pipe aperture decreases to the half of its initial aperture. On the other hand, when the normal contact force is tensile, the aperture of fluid channel can be obtained by:

(10)
where ͞ m is a dimensionless multiplier, and g refers to the normal distance between two particles.

For the fluid, it circulates into/out of the domain through the surrounding pipes (Fig. 3b). In each pipe, the flow rate q can be calculated using the Eq.8. During a given time increment Δt, the variation of fluid pressure in the domain can be identified by: The flow rate of each pipe is further updated in accordance with Eq.11. Furthermore, it is important to note that the fluid present within

0 0 0 n n n a F a F F = + 0 n F 0 a a mg = + ( ) f d d K p q t V V D = D -D å
the domain exerts pressure on the surfaces of the surrounding particles (Fig. 3b).

Improvement of pipe flow model

As mentioned in Eq.8, the seepage processes in the material depend strongly on the pipe aperture.

In the present work, the pipe is used to simulate the throat, with the pipe aperture representing the inner diameter of the throat. In natural state rock, the throat size is not uniform due to different geological sedimentation processes. However, the heterogeneous distribution of the initial pipe aperture is ignored in the literature. To address this issue, the present work proposes a new algorithm to take into account the heterogeneity of the pipe aperture.

Meanwhile, the core part of the pipe flow algorithm lies in the updating of fluid pressure and volume at each time step. The fluid flow through each pipe is established based on the cubic law (Eq.8), which has been confirmed by many studies. However, further discussion is required regarding the updating of the fluid pressure (Eq.11) within the domain. Currently, the pressure change is calculated based on the volume change of the fluid during each time increment Δt, assuming that the pressure change is totally induced by fluid compression. In practice, the pressure change is also related to the deformation of rock skeleton.

In order to take into account these previously mentioned phenomena, the previous pipe flow algorithm (Eq.11) will be improved in the present work. We begin with the non-uniform distribution of initial pipe aperture.

Inhomogeneous distribution of initial pipe aperture

Since the Weibull distribution function has been widely used to characterize the heterogeneity and anisotropy of geomaterials, it is adopted in the present work:

(12)
The probability density function of Weibull distribution is given by:

(13)
where x represents a random variable and refers to the initial aperture ɑ0 in the present work. η and β are the scale and shape parameters of the distribution, respectively. η controls the mean value of the random variable while β characterizes its distribution shape. When β ≤ 1, the density function is monotonically decreasing. For β > 1, the density function is unimodal with the mode ɑ0 = η{(β -1)/β} 1/β , and β reflects the dispersion degree of the random variable ɑ0. When β approaches infinity, the distribution of ɑ0 is close to that of the scale parameter. Conversely, smaller shape parameters indicate more dispersed distributions of the random variable ɑ0.

Previous experimental observations [START_REF] Feng | An improved geomechanical model for the prediction of fracture generation and distribution in brittle reservoirs[END_REF]) have shown that the aperture of microcrack vary from 0.5𝜇𝑚 to 12𝜇𝑚 in sandstones. Since η mainly controls the mean value of the fracture aperture, the value of η is taken as 1, corresponding to a mean fracture aperture of approximately 4 𝜇𝑚. In order to identify the parameter β, a series of preliminary study have been performed using three different β (β =2, 5 and 10). The statistics analysis of the number of pipes with different apertures is shown in Fig. 4. When β = 10, the aperture values vary within a very small range, predominantly concentrating around the mean value. When β = 5, the variation range of aperture values is about twice that of β = 10. Inversely, in the case of β = 2, the peak aperture value deviates from the mean value, and a significantly larger variation range compared to the other two cases.

Moreover, the spatial distributions of pipe apertures in the samples for these three values of β are ( )

1 x / F( x ) e b h - = - ( ) 1 x/ x f ( x ) e b b h b h h - - ae ö = ç ÷ è ø
shown in Fig. 5, where the line width represents the pipe aperture. When β =5, some thin lines can be found (Fig. 5b) and the overall variation of line width is slightly larger than that of β =10 (Fig. 5c). On the other hand, within the domain:

(14)
During each time step, the variation in pore pressure Δp is updated. According to the principle of action and reaction, it can be calculated based on the fluid deformation or by assessing the deformation of the rock skeleton. As a result, the variation in pore pressure Δp can be described as follows:

d f q t V V D = D + D å (15)
where 𝐾 # is the compression modulus of the rock skeleton. By combining Eq.14 and Eq.15, the volume change of the domainΔVd can be written as below:

(16)
Consequently, the variation in pore pressure Δp can be obtained by substituting Eq.16 into Eq.15:

(17)
Finally, the new pressure update formula is obtained by considering the deformation of rock skeleton.

The previous equation (Eq. 17) indicates a close relationship between the volume of the domain Vd and the pore pressure. In the literature, some researchers simplify the calculation by assuming the apparent volume of each domain to be the unit volume. However, this assumption is obviously unreasonable since the pore volume is much smaller than this value. To address this issue, previous studies (Chen et al. 2020) have used the difference between the volume of the model and the volume of all particles to estimate the domain's volume. They have also taken the mean value of these differences as the apparent volume of each domain. However, to consider the heterogeneity of the rock, a new approche is take in the present work: the volume of each domain is individually obtained through vector calculation based on the positions of the particles surrounding the domain.

Parameter identification and numerical validation

In this section, we focus on the identification of DEM parameters, including mechanical and hydraulic ones. The mesoscopic strength parameters of the contact model can be identified by using

f d f s d d V V p K K V V D D D = = f d s f K V qt K K D = D + å 1 1 1 d s f p qt V K K D = D ae ö + ç ÷ ç ÷ è ø å
a triaxial compression test while the hydraulic parameters will be determined via a steady-state seepage experiment. To calibrate these parameters, the laboratory tests performed by [START_REF] Xu | Permeability evolution of sandstone under short-term and long-term triaxial compression[END_REF] are used. These sandstone samples were taken from Yuncheng coal mine located at a depth of 960m. They are composed of gray fine sandstone with some dark red mudstone cementation. The volumetric density ranges from 2533 to 2591kg/m 3 and the permeability is approximately on the order of 12×10 -19 m 2 .

In the improved flow model presented in the previous section, the mesoscopic strength parameters do not appear in the seepage equations. Therefore, they do not have an impact on the seepage processes in the studied materials. As the emphasis of the present work is put on the influence of fracture properties on the seepage characteristics of rock mass, the mesoscopic strength parameters of the contact model are only identified by using one triaxial compression test. At the beginning, the sample without fractures is used to identify the mechanical parameters. After that, the improved pipe flow algorithm is implemented in PFC2D, and the hydraulic parameters of the intact sample are determined by using the steady-state seepage experiment.

Based on the empirical relations between meso-and macro-parameters in the particle discrete element method, the meso-parameters of the contact model are calibrated using trial-and-error method and given in Table 1. The numerical results and experimental data are compared in Fig 7a . A good agreement is obtained between the numerical and experimental results.

Subsequently, the hydraulic mesoscopic parameters are identified by the best fitting of the permeability of studied material. In the numerical simulation, the sample is free of mechanical loading.

Two water pressures of Pin = 1MPa and Pout = 0MPa are applied on the upper and lower boundaries of the sample, respectively. Based on the experimental works of [START_REF] Xu | Permeability evolution of sandstone under short-term and long-term triaxial compression[END_REF], an assembly of 100mm in width and 100mm in height should be adopted. To avoid the influence of boundary conditions and achieve a satisfactory estimation within a reasonable time, an assembly of 200mm in gradually increases. Eventually, the inflow and outflow converge towards the same value and a stable seepage sate is achieved. As a result, the permeability k of the sample can be calculated as follows:

(15

)
where Q is the steady flow rate, μ is the fluid viscosity, A and H are the width and height of the initial model, respectively. The seepage flow curves are shown in Fig. 7b. The computed steady-state seepage rate of the numerical sample is about 1.20×10 -9 m 3 /s, which corresponds to a permeability of 12×10 -19 m 2 . This permeability is basically consistent with that obtained by [START_REF] Xu | Permeability evolution of sandstone under short-term and long-term triaxial compression[END_REF]. The values of characteristic parameters are given in Table 1 and will be used to study the seepage characteristics of fractured samples, with a special attention to the influence of fracture aperture, fracture inclination and fracture length. According to the previous presented pipe flow model, the seepage process of rock mass depends strongly on the fracture properties. To characterize and comprehend the effects of non-penetrating fractures on the seepage behavior of rock mass, a series of parametric studies will be performed, using the samples fractured by a single non-penetrative embedded fracture, with varying fracture apertures 𝑎 " , fracture inclinations 𝜃 " and fracture lengths 𝐿 " .

To numerically reproduce the fracture, the contacts within a 1mm width surrounding the fracture geometry are selected to reproduce the microscopic properties of the fracture. Moreover, to facilitate the result comparison, the fracture center is positioned at the center of the sample (Fig. 8). The experimental observations in uniaxial/triaxial compression tests [START_REF] Feng | An improved geomechanical model for the prediction of fracture generation and distribution in brittle reservoirs[END_REF]) have shown that the length of fractures varies from several micrometers to a few centimeters. To get reliable numerical results and take into account the ratio of fracture size to sample size, three values of fracture length are adopted: Lf =6cm, 10cm and 14cm. In practice, the aperture of fracture generally varies within a large range, often being dozens or even tens thousands of times larger than the initial micro-cracks. However, for the computational efficiency reasons, the maximum fracture aperture is only set to be 100 times of that of the matrix in the present work. Three ratio of fracture aperture to its initial aperture (ɑf /ɑ0) of 10, 50 and 100 are adopted.

In geotechnical engineering, fractures are distributed randomly in the host rock, and the dip angle of fracture is an important issue in the seepage processes of rock mass. Five inclination angles of 0°, 30°, 60°, 90° and 120° are considered in the present work. All the numerical simulations are summarized in the Table 2. Visual presentations of some cases are given in Fig. 8. In practice, since the fracture is inclinedly distributed in the rock mass, we have taken the inclination angle of 𝛽 " =120°

as the reference angle while ɑf /ɑ0=50 and Lf =10cm are the reference aperture and length of fracture respectively. In order to quantify the influence of different fracture properties on the seepage characteristics of samples, eleven monitoring points are chosen in each fractured sample (Fig. 9). Among them, three points are located at the upper (called F1) extremity, middle (called F2) and lower extremity (called F3) of the fracture, respectively. The remaining points are placed within the non-fractured zone.

Using the fractured rock samples with different fracture properties, a series of seepage tests are simulated with the adoption of the same test scheme presented in the previous section (Fig. 6). The numerical results with and without fractures are firstly discussed to illustrate the role of fracture in the seepage processes of rock. After that, the influence of fracture properties will be discussed.

Effects of embedded fracture

The role of a non-penetrating embedded fracture in the seepage processes of rock mass is firstly investigated. Two calculations have been compared: with (called Ref) and without (called Null) fracture. The numerical results of two cases will be compared and analyzed in this section.

The curves of flow rate in the samples with and without fracture are given in Fig. 10. The general evolutionary trends of flow rates in two cases are the same: the inflow (𝑃 $% ) decreases progressively with time while the outflow (𝑃 &'( ) increases gradually, and finally these two flows tend to the same value. By comparing the obtained results with and without fracture, one observes that the flow rate obtained in the sample with fracture (Case Ref) is higher than that of sample without fracture (Case Null), indicating a larger permeability of the fractured sample (Table 2). Moreover, due to the presence of the fracture, more time is needed to reach the hydraulic equilibrium in the fractured sample.

The final distribution of pore pressure is shown in Fig. 11. It is observed that the distribution of pore pressure is significantly influenced by the presence of the embedded fracture. In the sample without fracture, the evolution of pore pressure is only observed in seepage direction, which is parallel to the vertical direction (y) of the sample. For a given height (i.e., y coordination), the same pore pressure is obtained. On the other hand, with respect to the fracture position, the pore pressure in the sample with fracture can be divided into two zones, called left and right zones. For a given y coordination, pore pressure in the left zone is higher than that of the right zone. These observations may be explained by the fact that the presence of fracture enhances the seepage processes of material, especially in the upper part where the higher pressure is imposed. The fracture acts as a conduit, accumulating water from the upper part and serving as a reservoir. As a result, the seepage process in the part below the facture is slowed down or decelerated.

Based on the obtained results, it can be seen that the embedded fracture has a non-negligible impact on the seepage behavior of rock. Therefore, it is crucial to further analyze the effects of fracture properties on the seepage characteristics of rock masses. By studying these effects, a more comprehensive understanding of how fracture properties influence seepage behavior of rock mass can be obtained.

Effects of fracture aperture

The influence of fracture aperture on the seepage characteristics of rock mass is discussed in this section. In the numerical simulations, the fracture length is 10 cm and the fracture dip angle is 120° while three fracture apertures ɑf /ɑ0 of 10, 50 and 100 are used.

As shown in Fig. 12, the inflow and outflow rates obtained in three samples are firstly compared.

The similar trend is observed in all cases. At the beginning, the inflow decreases progressively with time while the outflow gradually increases. After that, the inflow and outflow rates tend to converge to the same value. Before reaching steady state, the inflow rate increases with the fracture aperture while the outflow decreases with increasing fracture aperture. For instance, the sample with ɑf /ɑ0 =100 exhibits the highest inflow and the lowest outflow, while the sample with ɑf /ɑ0 =10 has the lowest inflow and highest outflow. Finally, at the steady-state phase, three flow rates of 1.4221×10 - 9 m 3 /s、1.4247×10 -9 m 3 /s、1.4253×10 -9 m 3 /s are obtained in the samples with ɑf /ɑ0 =10, 50 and 100, respectively. Although slight differences are observed among these three cases, the effect of fracture aperture on the seepage rate is still noticeable: the flow rate increases with the increase of fracture aperture and the most rapid hydraulic equilibrium is achieved in the case with the smallest af/a0, i.e.

af/a0 = 10.

Moreover, the permeability of the fractured specimens is compared and analyzed. It is observed that the permeability of the fractured rock sample increases with increasing fracture aperture: k = 14.220×10 -19 m 2 with af/a0 = 10, k = 14.247×10 -19 m 2 with af/a0 = 50 and k = 14.253×10 -19 m 2 with af/a0 = 100. This observation is coherent with the evolution of flow rate: the larger the flow rate, the greater the permeability. Similar results were also obtained in the experimental tests performed on the samples with penetrating fractures [START_REF] Qian | Experimental evidence of scale-dependent hydraulic conductivity for fully developed turbulent flow in a single fracture[END_REF]). However, as the fracture is non-penetrating and completely embedded in the sample, the influence of fracture aperture on the permeability of rock mass is significantly limited compared to that of penetrating fractures. It is important to note that the numerical results and mesoscopic analysis in the present work is capable of providing an estimation of the effects of embedded fractures on the seepage characteristic of rock mass. This contributes to enhancing our understanding on the role of non-penetrating fracture on the seepage behavior of rock mass. The distributions of pore pressure in the samples are shown in Fig. 13. Overall, the pore pressure distributions exhibits striking similarities among three cases: two seepage areas, namely left and right zones, can be generally observed (Fig. 13 a-c). For a given y value, the pressure in the right zone closely approximates that of the left one in the upper part of the sample. However, in the central and lower parts of the sample, the pressure obtained in the left zone is slightly higher than that of right zone. Moreover, the evolutions of pore pressure at different measured points are given in Fig. 13d. It seems that the increase of fracture aperture has a minimal effect on the pore pressure distribution in the sample.

To quantitatively assess the influence level of fracture aperture on the evolution of pore pressure within the fracture, the pressure evolutions at monitoring points F1, F2 and F3 are given in Fig. 14. It can be observed that the rate of pore pressure increase within the fracture decreases with increasing fracture aperture. The smallest increasement of pressure within the fracture is observed in the case with the biggest af/a0, i.e af/a0 = 100. Moreover, the hydraulic equilibrium is delayed by the increase of fracture aperture: the final equilibrium of pore pressure is reached more slowly in the sample with larger fracture apertures (af/a0 = 50, 100) compared to the samples with a small fracture aperture (af/a0 = 10). After that, the pressures obtained during the steady-state phase are given in Fig. 15. In Fig. 15a, the results of three monitoring points (F1, F2, and F3), located within the fracture, are compared. In general, for a given fracture aperture, the pressure gradually decreases from the top to the bottom of the fracture (from F1 to F3), and the pressure difference between the top and bottom of the fracture (between F1 and F3) decreases with increasing fracture aperture. For instance, when af/a0 = 10, 50 and 100, the pressure difference between F1 and F3 points is 8kPa, 4kPa, and 3kPa, respectively. Moreover, for a given case, the pressure observed at the upper extremity of fracture (F1) decreases with increasing fracture aperture, while the pressure obtained at the lower extremity of fracture (F3) increases with increasing fracture aperture. The pressure in the center of the fracture (F2) remains nearly constant in different cases.

As mentioned earlier in Fig. 13d, the influence of fracture aperture on pore pressure evolution in the non-fractured zone is limited. To evaluate the influence of fracture aperture on pore pressure in the non-fractured zone of the sample, the relative pressures at the observation points within the nonfractured zone are compared in Figs. 15b-d. For a given point, the relative pressure is obtained by subtracting an identical reference value from the obtained results. As the reference values are different for the different monitoring points, the obtained pressures only indicate the effects of fracture aperture on the pressure evolution at the given measurement point. At two points S1 and S8, located in the farfield of the fracture, the relative pressure remains largely undisturbed by the variation in fracture aperture. On the other hand, at monitoring points S4, S6 and S7, the relative pressure increases with increasing fracture aperture, while the opposite trend is observed at monitoring points S2, S3 and S5.

It appears that a dividing line exists in the sample: with the increase of fracture aperture, pore pressure decreases in the zone above this line while pore pressure increases in the zone below this line. For instance, in the specimen with a fracture inclination of 120°, the dividing line may align with the direction that is centrosymmetric to the fracture direction and pass through the center of the specimen.

The similar dividing line was also identified in the specimen with a fracture inclination angle of 30°.

Effects of fracture dip angle

The numerical results of the samples with different fracture dip angles are compared and discussed in this section. In the numerical simulations, the fracture length is Lf = 0.1m and the fracture aperture is ɑf /ɑ0 =50 while the fracture dip angle βf varies from 0° to 120°.

The evolutions of flow rates in the samples with different βf are given in Fig. 16. When βf increases from 0° to 90°, both the outflow and inflow rates progressively increase. Since the influence of the fracture inclination angle on the seepage characteristics is primarily related to the angle between the fracture inclination and the seepage direction, the numerical results of the samples with the fracture dip angle ranging from 0° to 90° should be equivalent to those with dip angle decreasing from 180° to 90°. In the present work, it is also noticed that the evolutions of flow rates in the sample with βf = 120° are nearly identical to those with βf = 60°.

The permeability of fractured rock mass is also affected by βf . When βf increases from 0° to 120°

(with a constant increment of 30°), the permeability of rock varies as follows: 11.81×10 -19 m 2 , 12.92×10 -19 m 2 , 14.31×10 -19 m 2 , 15.52×10 -19 m 2 , 14.25×10 -19 m 2 , respectively. It can be noticed that with the increase of the fracture inclination angle (from 0° to 90°), the permeability progressively increases. The two samples with βf = 120° and 60° exhibit nearly the same permeability values although there is still a slight difference between them. This difference may be related to the heterogeneity of the initial sample. Once again, the consistent results obtained from the two samples with βf = 120° and 60° confirm the reliability of the proposed algorithm.

In the steady-state phase, the distributions of pore pressure in the samples with different βf are given in Fig.17 significantly affected by βf. Notably, in two samples with fracture tip angles of 0° and 90°, a symmetrical distribution of pressure with a horizontal or vertical line shape is observed. Moreover, for the samples with an inclination angle between 0° and 90°, the pressure distribution exhibits a stepped shape.

To quantitatively assess the influence of βf on the pressure evolution, the pore pressures at the monitoring points are compared in Fig. 18. Firstly, the results of three points located within the fracture (F1, F2 and F3) are analyzed (Fig. 18a). In all cases, the pressure generally decreases from the top (F1)

to bottom (F3) of fracture, except for the case with a horizontal fracture where the pressure remains almost the same at three observation points within the fracture. As the inclination angle of the fracture increases from 30° to 90°, the pressure difference between the upper and lower extremities of the fracture gradually increases.

Moreover, the pore pressures in the non-fractured zone of the samples are shown in Figs. 18b-d.

To investigate the influence of fracture dip angle, the numerical results of the sample without fracture (called Null) are also added in the comparison. It is important to note that the positions of the fracture extremities vary due to the changing fracture inclination. In the upper region (S1-S3 points) of the specimen, the pore pressure at the S2 point, located in the near-field of the upper extremity of the fracture, is strongly disturbed by βf : the pressure decreases with increasing βf when 0° < βf < 90°.

Conversely, for the other two points (S1 and S3), located in the far-filed of the upper extremity of the fracture, pressures are less disturbed by the variation of βf. For instance, when βf increases from 0° to 90°, the pressure variations obtained at S1 and S3 are only approximately 30kPa while the pressure variation at S2 is about 100kPa. Comparing the pressures obtained in samples with different βf, it can be observed that the pressure of point S2 is the smallest in the sample with a vertical fracture (βf = 90°) among all the monitoring results of point S2. This observation is due to the fact that S2 is the closest point to the upper extremity of the fracture in the sample with βf =90°. Moreover, a similar trend can be found in the lower zone (S6-S8) of the specimen (Fig. 18d). When βf increases from 0° to 90°, a variation of pressure about 30kPa is observed at the two points (S6, S8) while the pressure variation attains 70kPa at S7, which is the point closest to the lower extremity of fracture compared to S6 and S8. Meanwhile, at S7, located in the near field of the fracture, the pressure progressively increases with increasing βf with the range of 0° < βf < 90°, and the maximum pressure is obtained in the sample In contrast, in the central part of the sample, the pressures at the monitoring points S4 and S5 are more disturbed by βf compared to F2, which is located in the center of the sample. These observations confirm that βf plays a significant role in pore pressure evolution within the fractured sample.

Additionally, it is noteworthy that in all the fractured samples, the pore pressure obtained at F2 is smaller than that of the non-fractured sample while the pressures at different positions in the specimens with a horizontal fracture are almost the same as those of the sample without fracture, except for the monitoring point F2.

Effects of fracture length

In this section, the effects of fracture length on the seepage characteristics of rock samples are further investigated. A series of seepage tests are realized with fracture length Lf = 0.06m, 0.1m, and 0.14m while the other fracture properties remain the same.

As shown in Fig. 19, the seepage flow rate of rock sample is significantly affected by the fracture length. The permeability values of the samples with the fracture lengths of 0.06m, 0.10m and 0.14m are measured to be 12.47×10 -19 m 2 、14.25×10 -19 m 2 and 16.18×10 -19 m 2 , respectively. It is worth nothing that the permeability of the sample increases progressively with increasing fracture length.

The final distributions of pore pressure in the specimens are given in Fig. 20. Since the inclination angles of fracture in three samples are the same, a stepped shape is generally observed in the pore pressure distribution. Moreover, the step size seems to increase with increasing Lf. For instance, in the case where Lf = 0.06m, the influence of the fracture on the pore pressure distribution is primarily restricted in the central part of the sample. When the fracture length reaches 0.14m, the disturbance in pore pressure is observed in a larger portion of the sample experiences due to the presence of the fracture.

The pressures of all monitoring points are given in Fig. 21. It is observed that the pressure within the fractures increases as the length of the embedded fracture increases (Fig. 21a). The average pressure in the fracture with Lf = 0.14m is approximately 560kPa while the average pressures in the fractures with the lengths of 0.06m and 0.1m are around 500 and 510kPa. The pore pressures obtained at the observation points in the non-fractured zone are given in Fig. 21b-d. The pressures at points S1 and S8, located in the far-filed of the fractures, remain largely unaffected by variation in Lf. However, at points S2, S3 and S5, the pressures decrease progressively with an increase in fracture length.

Moreover, at a points S2, S3, and S5, the pressures are lower compared to those of the sample without fracture. On the other hand, at points S4, S6 and S7, the pressures increase with increasing of fracture length and are higher than those of the sample without fracture. Once again, a dividing line may exist in the sample: with the increase of fracture length, pore pressure decreases in the zone above this line while it increases in the zone below this line. In the specimen with a fracture inclination of 120°, the divided line may be in the direction that is centrosymmetric to the fracture inclination and passes through the center of the specimen.

In the center of the specimen (F2), the pore pressure in the sample with the fracture length of 0.14m is about 30kPa higher than that of the non-fractured specimen while the pore pressure of nonfractured specimen is approximately 30kPa higher than that of the samples with fracture length of 0.1m and 0.06m. This observation may be due to the fact that the seepage mechanism of the singlefracture specimen changes when the fracture length exceeds a certain value. consistent with previous research works [START_REF] Qian | Experimental evidence of scale-dependent hydraulic conductivity for fully developed turbulent flow in a single fracture[END_REF][START_REF] Singh | Effect of sample size on the fluid flow through a single fractured granitoid[END_REF][START_REF] Wu | Investigations on the seepage characteristics of fractured sandstone based on NMR real-time imaging[END_REF]).

However, while previous studies often attributed a similar level of influence to fracture aperture, dip angle and length on the seepage process of rock, the present study reveals two levels of influence on the seepage characteristics of rock. Specifically, the fracture length and dip angle have an important impact on the seepage characteristic of rock mass while the influence of embedded fracture aperture is very relatively limited. This difference may be due to the nature of the embedded fracture, which differs from the penetrating fracture.

On the other hand, the pore pressure distribution in the rock mass also depends on fracture properties. The effects of fractures on pore pressure distributions can be explained by the hysteresis of pore pressure growth within the fractures. In this regard, the distributions of pore pressure at different seepage stages in the sample with a fracture inclination angle of βf = 120° is discussed (Fig. 22). At the beginning (Fig. 22a), the pore pressure in the upper part of the sample increases progressively, with similar pressure growth observed for a given height. After that, as fluid continues to penetrate the sample, the pressure front reaches the upper extremity of fracture (Fig. 22b). However, due to the long and narrow nature of the fracture, it acts as a cavity that accumulates fluid, leading to a gradual increase in internal pressure as the pressure is transmitted from the nearby pore zones to the fracture. As a result, for a given height, the pressures in the near field of the fracture cannot increase synchronously with the pressure in the far field of the fracture.

Meanwhile, with respect to the fracture position, the upper part of the specimen is divided into two distinct pressure transfer zone (Fig. 22c): a left seepage zone (referred to as zone Ⅰ) and a right seepage zone (referred to as zone Ⅱ). Given that the inclination angle of the fracture is 120°, the cross-sectional area of zone Ⅰ is considerably larger than that of zone Ⅱ. Consequently, a greater amount of pore water is transmitted to the lower part of the sample in zone Ⅰ, while only a small amount of pore water is transmitted to the lower part of zone II. As a result, the pressure growth in the left and right seepage zones of the sample is not synchronized, and the pore pressure in the zone I increases at a slower rate compared to zone II. Consequently, in the samples with βf = 30° and 60°, higher pressure is obtained in the right zone with respect to the left zone while the pressure is uniform in both the left and right parts of the sample when βf = 90° . (1) The permeability of fractured rock mass with non-penetrating embedded fracture depends strongly on the dip angle and length of non-penetrating fracture and is slightly impacted by the fracture aperture. In general, permeability increases with increasing fracture aperture and length, and decreases with increasing angle between the embedded fracture and the seepage direction.

(2) The effects of non-penetrating embedded fracture on the rock permeability are different from those of through-penetrating one. Firstly, the penetrating fractures can induce an increase of several orders of magnitude of rock permeability while the permeability of rock containing non-penetrating fractures is almost of the same order as that of non-fractured rock. Secondly, in the rock mass with non-penetrating fractures, the influence of the aperture of embedded fracture on the permeability of rock is less important than that of penetrating one.

(3) The pore pressure is slightly disturbed by the variation of fracture aperture and depends strongly on the dip angle and length of fracture. A dividing line seems to exist in the sample:

with the increase of fracture length/aperture, pore pressure decreases in the zone above this line while pore pressure increases in the zone below it.

(4) The average pore pressures inside the fractures are relatively close when the fracture length is smaller than 0.14m. It seems that the seepage mechanism in the sample fractured by a nonpenetrative embedded fracture changes once the ratio between the fracture length and the sample dimension is bigger than 0.7. Further investigation is needed to understand this phenomenon.

Overall, the improved pipe fluid algorithm provides valuable insights into the influence of nonpenetrating embedded fractures on seepage characteristics. The numerical analysis and discussions can help to get a better understanding of rock mass behavior and provide a practical guidance in engineering and geotechnical applications.

  Fig.1 Force-displacement law of parallel bond model

Fig

  Fig.2 The fluid network model

  . The pipe aperture ɑ can be obtained by the following equation:

where

  Kf is the fluid bulk modulus, Vd is the domain volume, ΔVd refers to the volume change of the domain during the time step Δt. ƩqΔt represents the cumulative sum of fluid volume changes within the studied domain. Finally, by using Eq.11, the fluid pressure in each domain can be updated based on the calculated variation of fluid pressure. This update considers the fluid volume change within the domain during the time step Δt. Additionally, this update allows for the determination of a new pressure gradient between the two ends of each pipe associated with the domain. The pressure gradient is crucial for driving fluid flow through the pipes.

  Fig. 5a (β =2) clearly displays difference in line width: the think lines in lawn green are randomly distributed while the medium lines in simple green dominate and thin lines in dark green are also present. The heterogeneity of pipe aperture is satisfactorily reproduced with β = 2. Therefore, η =1 and β = 2 are adopted in the present work to describe the non-uniform distribution of the initial fracture aperture ɑ0.

Fig. 4

 4 Fig. 4 Statistics of the number of pipes with different apertures

Fig. 6

 6 Fig.6 Steady-state seepage model

  Fig.8 Numerical samples with different fracture properties

Fig. 10

 10 Fig.10 Evolutions of flow rate in the samples with fracture and without fracture

Fig. 12

 12 Fig.12 Flow rate curves of samples with different fracture apertures

  Fig.13 Pressure distributions in the samples with different fracture apertures

  29

Fig. 16

 16 Fig.16 Evolutions of flow rate in the samples with different fracture dip angles

  Fig.17 Pore pressure distribution of samples with different fracture dip angles

  Fig.18 Pore pressure of monitoring points in specimens with different fracture dip angles

Fig. 19

 19 Fig.19 Flow rate curves of samples with different fracture lengths

  Fig.22 Pore pressure distribution in the reference sample at different seepage stages

Table 1

 1 Mechanical and hydraulic parameters used in the numerical simulation

	Property	Value
	Bond effective modulus E * /GPa	14.5
	Effective modulus E /GPa	14
	Normal-to-shear stiffness ratio k *	2.5
	Friction coefficient μ	1.5
	Bond gap g/m	1e-5
	Radius multiplier λ	0.8
	Moment contribution factor β	0.5
	Tensile strength σc /MPa	51
	Cohesion c /MPa	51
	Friction angle φ /°	50
	Particle bulk modulus, Ks / GPa	20
	Fluid viscosity μ / Pa s	1e-3
	Fluid calculation timestep Δt / s	1e-3
	Initial hydraulic aperture ɑ0 /m	4.4e-7
	Fluid bulk modulus, Kf / GPa	2

modelling the sandstone samples with a single embedded fracture

Table 2 .

 2 Parametric study cases

	Influence	Cases	Fracture aperture 𝑎 ! 𝑎 " ⁄	Fracture angle 𝛽 ! (°)	Fracture length 𝐿 ! (m)	Permeability k (𝑚 # )
	Intact sample	Null	0	0	0	12×10 -19
	Reference	Ref	50	120	0.1	14.25×10 -19
	Fracture aperture	1 2	10 100	120	0.1	14.22×10 -19 14.25×10 -19
		3		0		11.81×10 -19
	Fracture angle	4 5	50	30 60	0.1	12.92×10 -19 14.31×10 -19
		6		90		15.52×10 -19
	Fracture length	7 8	50	120	0.06 0.14	12.47×10 -19 16.18×10 -19

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (Grants Nos. 41831278), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX22_0615) and the National Key R&D Program of China (2018YFC1508501).

Data Availability statement

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

ORCID

Xiao Chen: https://orcid.org/0000-0003-3250-7319

Chong Shi: https://orcid.org /0000-0003-1386-0651 Yun Jia: https://orcid.org/0000-0002-2844-8777

on the seepage characteristics of fractures [START_REF] Wu | Investigation of the effective hydro-mechanical properties of soil-rock mixtures using the multiscale numerical manifold model Computers and[END_REF][START_REF] Wu | Numerical manifold computational homogenization for hydro-dynamic analysis of discontinuous heterogeneous porous media[END_REF], most of these studies treated the fracture and the solid phase as two independent entities rather than considering them as a unified system. In contrast, the present work considers the fracture and rock matrix as an integrated fractured material, allowing us to accurately capture the effects of fracture properties on the seepage characteristics of rocks using an improved fluid algorithm. The aperture, dip angle and length of the non-penetrating embedded fracture have a significant impact on the permeability and pressure distribution of the rock mass.

Generally, the permeability of rock mass progressively increases with increasing the aperture and length of the non-penetrating fracture. Conversely, the rock permeability decreases as the angle between the embedded fracture and the seepage direction increases (Table 2). These observations are

Author contributions

Xiao Chen: Methodology, numerical calculation, writing-original draft, writing-review & editing, funding acquisition. Chong Shi: Software, methodology, funding acquisition, supervision.

Yiping Zhang: Writing-review & editing, software. Yun Jia: Supervision, writing-review & editing.

Huaining Ruan: Funding acquisition, project administration.