The equivalent discrete fracture networks based on the correlation index in highly fractured rock masses

Engineering Geology 260 (2019) 105228

Contents lists available at ScienceDirect

Engineering Geology journal homepage: www.elsevier.com/locate/enggeo

The equivalent discrete fracture networks based on the correlation index in highly fractured rock masses

T

⁎

Guowei Maa,b, Tuo Lib, , Yang Wangb, Yun Chenb a b

School of Civil and Transportation Engineering, Hebei University of Technology, Tianjin 300401, China School of Civil, Environmental and Mining Engineering, The University of Western Australia, Perth, WA 6009, Australia

A R T I C LE I N FO

A B S T R A C T

Keywords: Highly fractured rock masses Correlation index Permeability similarity Equivalent permeability factor Density-reduced models

In the numerical simulations of highly fractured geological formations, discrete approaches are considerably promising and adequate to describe fluid flow in detail. However, the computational complexity increases dramatically with a greater number of fractures. This becomes the primary limitation for field-scale applications. In this study, a correlation index is for the first time introduced to evaluate the significance of individual fractures, and an equivalent model is proposed to mimic the original domain with a density-reduced one. By an equivalent permeability factor, the suggested model simplifies computational complexity, but compromises result precision to minor extent. This approach is validated in typical discrete fracture networks generated with stochastic fractal models. Effects of fracture geometry are discussed based on various distribution patterns. This method improves mesh quality when dealing with a fracture-matrix domain. It is also capable of optimizing reservoir design through fast and accurate estimations of gas productivity under different boundary conditions.

1. Introduction Rock masses are lithologically heterogeneous and anisotropic. They contain discontinuities of various scales and types, including cracks, fissures, joints, faults, and fault zones (generally referred to as fractures). Fractures embedded in rock masses consist of the principal conduits for a fluid flow (Liu et al., 2016; Pan et al., 2010). Deriving accurate and efficient simulations of fluid flow is the prerequisite of multi-physics coupling applications, such as hydraulic fracturing and pollutant migration (Ma et al., 2019a; Zhou et al., 2018). The most widely accepted numerical approaches to model fracture networks are conceptually categorized into continuum and discrete models (Ko et al., 2015; Larsson et al., 2013; Lee et al., 1995; Maryška et al., 2005; Zhou et al., 2018). Continuum models is preferred in handling large-scale problems, but their assumption of homogenization overlooks the geometric and hydraulic properties of individual fractures (Liu et al., 2018; Reeves et al., 2013). In naturally fractured reservoirs (NFRs), discontinuities are highly random and heterogeneous in size, location, orientation, and aperture. The inherent assumption of uniformity in the continuum models generates significant discrepancies between numerical simulations and reality. In contrast, the discrete models may appear more logical since they reduce the abstractions and explicitly incorporate each fracture in the domain (Chen et al., 2019b; Ma et al., 2017, 2019b). According to statistics from site investigations, fractures ⁎

are commonly approximated by line segments in two dimensions (Hamdia et al., 2017) or planar disks in three dimensions (Chen et al., 2018a). Although explicit representations of fractures retain result precision to some extent, they lead to a significant computational burden, especially in a field-scale problem. Three major factors contribute to the computational complexity of discrete methods: (1) the significant quantity of discontinuities, (2) the uncertainty of data from site investigation, and (3) the constraints of the discretization in a fractured porous medium. The number of fractures is the primary contributor to the computational complexity in discrete methods. The density reduction of fractures has been widely discussed in discrete fracture networks (DFN) models. They neglect flow in rock matrix and only concern the contribution of fracture networks as rocks have high contrast in permeability. Conventional methods reduce the number of fractures by neglecting isolated fractures and trimming the nonflowing dead ends (Bidgoli et al., 2013; Parashar and Reeves, 2012). However, the effectiveness of this technique decreases dramatically in highly fractured rock masses, where fracture networks are well connected. When the computational capacity is incapable of handling numerous discontinuities explicitly, compromising between accuracy and efficiency is necessary. Ren et al. (2017) proposed the equivalent discrete fracture networks (E-DFN) method to reduce the simulation complexity by sacrificing some extent of result accuracy. This method improves

Corresponding author. E-mail address: [email protected] (T. Li).

https://doi.org/10.1016/j.enggeo.2019.105228 Received 15 August 2018; Received in revised form 26 May 2019; Accepted 8 July 2019 Available online 09 July 2019 0013-7952/ © 2019 Elsevier B.V. All rights reserved.

Engineering Geology 260 (2019) 105228

G. Ma, et al.

Abbreviations DFN EDFM E-DFN LSD NFR REV TPFA

E Young's modulus E0, …, E90 Equivalent domain with corresponding fracture-density reduction Fi, Ri A random number in the range of [0, 1] K Fisher constant KIC Stress intensity factor Ln Domain size Mr Enveloped area of the permeability curve of density-reduced domain with reduction rate r N Normalized intersection number of fracture Pc Density probability of fracture Ps Probability of fracture centers in the corresponding subdomain Pc Normalized density probability of fracture R(l) Cumulative probability function of fracture length R0, …, R90 Density-reduced domain with corresponding fracturedensity reduction S Correlation index of fracture T Total number of subdomains in a parent domain in a multiplicative cascade process α Normalization factor in a power law distribution β Significance factor of fracture density γ Significance factor of fracture aperture μb Mean value in a log-normal distribution of fracture aperture σb Standard deviation in a log-normal distribution of fracture aperture εr Variance in similarity between the equivalent and original permeability curves ϕi Individual fracture orientation ϕ Mean value of fracture orientation in a Fisher distribution

Discrete fracture networks Embedded discrete-fracture modeling Equivalent discrete fracture networks Limit state design Naturally fractured reservoir Representative elementary volume Two-point flux approximation

Symbols a b b krθ

krθ l li lmax, lmin n(l) p q sr v ωr Dq

Length exponent in a power law distribution of fracture length Fracture aperture Normalized hydraulic aperture of fracture Average permeability of density-reduced domain with reduction rate r and testing direction θ Normalized average permeability of density-reduced domain Fracture length Individual fracture length Upper and lower bounds of fracture length Probability density function of fracture length Percolation coefficient in a power law distribution of fracture length Dimension of order in a fractal model of fracture centers Length ratio between the side of a parent domain and its subdomain Poisson's ratio Equivalent permeability factor of domain with density reduction rate r Fractal dimension in a fractal model of fracture centers

permeability of fracture networks and included the effects of flow tortuosity. However, these models depend on the existence of the representative elementary volume (REV). When the simulation domain shows dramatic heterogeneity, the results obtained based on the assumption of REV may be unreliable. The constraints of discretization add another computational burden to discrete models. When modeling a fractured porous medium, most discrete methods rely on the high quality of meshes (Wang et al., 2017). Because of the constraint of conformity, the number of meshed elements and the quality of discretization are significantly influenced by the ill-conditioned intersections. The improvement in mesh quality leads to higher computational costs at the preprocessing stage and larger degrees of freedom at the calculation stage. Field-scale applications, such as oil extraction (Chen et al., 2019a) and geothermal exploitation (Chen et al., 2018b; Chen et al., 2019b), commonly involve hundreds of fractures, and hence up to a million computing nodes. This study proposes a discontinuity equivalence approach that extends the density-reduced technique to more general engineering practices. It for the first time defines a correlation index based on the geometric configuration of fracture networks. This index evaluates the significance of individual fractures to the domain transmissibility quantitatively. The overall flow rate in a fractured medium is dominated by the connectivity between boundaries and the hydraulic aperture of individual fractures (Hamm et al., 2007; Le Borgne et al., 2006). In literature (Berkowitz, 1995; Huang et al., 2017; Ren et al., 2015), the connectivity in random fracture geometries is measured by the average number of intersections per fracture. The contribution of individual fractures to the domain connectivity is positively correlated with the fracture length (Darcel et al., 2003; Renshaw, 1996; Xu et al., 2007). Another contributing factor that determines the fracture

computational efficiency through the reduction of fracture density. Nevertheless, E-DFN regards each fracture as an equally important component in the fracture networks. Fractures over different sectors are randomly disregarded and eliminated in terms of location, length, and aperture without evaluation. Moreover, the fracture networks generated in E-DFN are hardly realistic as they are identical in length and uniformly distributed across the entire domain. In NFRs, the configurations of fracture networks commonly follow fractal models, and the fluid flow in the domain is dominated by a group of critical conduits. When the fracture distributions become more realistic, the extent of density reduction by E-DFN is limited accompanied by a dramatic increase in errors and uncertainty. The uncertainty of fracture data deteriorates computational complexity. It is due to the fact that a practical means to use discrete models is to generate a series of stochastic realizations based on the data collected and conducting statistical analysis (Chesnaux et al., 2009; Earon and Olofsson, 2018; Hamdia et al., 2018). Vu-Bac et al. (2016) developed a Matlab toolbox for uncertainty and sensitivity analysis in computational expensive models. However, the usual practice in Monte Carlo simulations is assumed that each fracture is generated randomly based on statistics. In site investigations, the certainty of individual fractures varies from each other, which is related to their dimensions. We tend to have better confidence in the data of faults and joints compared with information related to fissures and cracks. Another possible approach to predict the permeability of stochastic fracture systems is through fractal characteristics in the fracture distribution of rock masses (Lu et al., 2016; Yang et al., 2017; Zhao et al., 2009). Yu et al. (2005) developed a two-dimensional fractal model that can predict the transport properties of saturated or unsaturated porous media. Liu et al. (2015) extended the fractal models to derive the equivalent 2

Engineering Geology 260 (2019) 105228

G. Ma, et al.

Fig. 1. Domains with different extent of fracture reduction: (a) R0; (b) R50; (c) R70; (d) R90.

restored from the density-reduced domain with an equivalent permeability factor. This technique enhances the computational efficiency of simulation and retains most discontinuous properties in fractured rock masses. It is not only applicable in complex DFN models but also preferable when preprocessing a fracture-matrix system. As the precision requirements may vary at different stages of a project, the proposed approach provide the quantitative basis for constructing the simulation model with various error tolerances. In this paper, we introduce the generation of the stochastic fractured domain with extreme heterogeneity and anisotropy. Then, a densityreduction process is introduced to a typical DFN model based on the correlation index. Multiple cases are tested and discussed to validate the improved equivalent model according to different fracture configurations. This model also demonstrates its applicability to practical problems, such as the simulations of fracture-matrix systems and reservoir-design optimization.

Table 1 Input parameters for typical DFN domain. Parameter

Symbols

Values

Units

Domain size Minimum fracture length Maximum fracture length Length exponent Percolation coefficient Constant aperture Fragmentation ratio Order of dimension for centers Fractal dimension of centers

Ln lmin lmax a p b sr q Dq

10 0.3 15 2.5 80 3 2 2 1.8

m m m − − mm − − −

Table 2 Input parameters for joint sets. joint set

1 2

Mean value (ϕ )

Fisher constant (K)

[deg]

[−]

45 135

25.0 35.0

2. Generation of fracture networks As reported in the literature, the fractures in two-dimensional rock masses are regarded as fractal-like tree networks. Hence, the cumulative size distribution of fracture length follows a power law distribution (Chen et al., 2018; Huang et al., 2017; Reeves et al., 2013), which yields

transmissibility is its hydraulic aperture (Bear, 1988; Li et al., 2016; Yin et al., 2017). The importance of individual fractures shall also be weighed by the effects of physical location, especially when fracture clustering exists. Based on this correlation index, the density of fracture networks is reduced to generate a DFN skeleton consisting of a few selected fractures. The hydraulic properties of the original domain are

n (l) = αl−a

l ∈ [lmin lmax ]

(1)

where n(l)[−] is the probability density for fractures with a length of l [L]; lmax[L] and lmin[L] are the upper and lower bounds of fracture length detected; a[−] is the length exponent; and α[−] is the 3

Engineering Geology 260 (2019) 105228

G. Ma, et al.

normalization factor. Accordingly, the cumulative probability R(l) [−] can be expressed as

R (l) =

∫l

l

n (l) dl =

min

α −a + 1 (l−a + 1 − lmin ) −a + 1

(2)

Hence, the normalization factor and the length of individual fractures li[L] can be derived as below:

α=

−a + 1 −a + 1 −a + 1 lmax − lmin

(3) 1 −a + 1

−a + 1 −a + 1 ⎤ ⎞Ri + lmin li = ⎡⎛ ⎣⎝ α ⎠ ⎦

(4)

where Ri[−] is a random number in the range of [0, 1]. The fault numbers in the simulation are roughly similar despite very different system sizes (Bour and Davy, 1997). However, the mass density of fractures in the domain is adjustable by the percolation coefficient p [−], which yields

p=

∫l

lmax

min

min(Ln2 l 2) n (l) dl Ln2

(5)

where Ln[L] is the domain size. The length exponent and percolation coefficient determine the system connectivity. The generation processes of fracture lengths by a fractal model have been introduced and discussed comprehensively in previous studies (Darcel et al., 2003; de Dreuzy et al., 2001; Liu et al., 2015). The aperture also significantly impacts the hydraulic properties of the fracture. Three types of distributions are discussed. The fracture aperture b[L] is assumed to be constant to highlight the effect of fracture length. Then the log-normal distribution (de Dreuzy et al., 2001; Pan et al., 2010) is applied through the typical probability density function as below:

Fig. 2. Domain setups in determination of directional permeability curves: (a) rotational simulation domains with 15∘ angle per step (dotted square); (b) boundary conditions of simulation domain.

Fig. 3. Similarity in permeability curves during fracture density reduction: (a) permeability curves of density-reduced domains; shape similarity in the equivalent domain of (b) E50, (c) E70, (d) E90. 4

Engineering Geology 260 (2019) 105228

G. Ma, et al.

Fig. 4. Effects of fracture aperture distribution: (a) constant; (b) length-related; (c) log-normal with μb = 5mm and σb = 1mm; (d) log-normal with μb = 5mm and σb = 2.5mm; (e) errors in permeability curves.

f (b, μb , σb) =

1 1 exp ⎧− 2 [ln (b) − μb ]2 ⎫ ⎨ ⎬ 2π σb ⎩ 2σb ⎭

b= (6)

KIC (1 − v 2) l E π /8 1

(7)

where KIC ⎡M·L− 2 ·T −2⎤ is the stress intensity factor and v[−] and E ⎣ ⎦ [M · L−1 · T−2] are Poisson's ratio and Young's modulus, respectively. Two types of distributions are adopted in the generation of fracture centers, namely, multiplicative cascade (Darcel et al., 2003; de Dreuzy et al., 2004; Meakin, 1991; Schertzer and Lovejoy, 1987) and Poisson

where μb[m] and σb[m] are the first and second moments of the lognormal distribution. Fracture apertures also depend on their length when subjected to in-situ stress (Baghbanan and Jing, 2007; Olson, 2003; Renshaw and Park, 1997; Zhang et al., 2017). In this case, the length-related aperture is expressed as 5

Engineering Geology 260 (2019) 105228

G. Ma, et al.

Fig. 5. Effects of fracture location distribution: (a) Poisson process; fractal models with (b) Dq = 1.95, (c) Dq = 1.85, (d) Dq = 1.75; (e) errors in permeability curves.

Fisher distribution is adopted in this study. The orientation of individual fracture ϕi[deg] can be generated through the equation as below:

(Liu et al., 2018; Priest, 1993; Ren et al., 2017) processes. The former process iteratively divides a domain into subdomains and assigns a probability Ps to each of them. For every parent domain, such probabilities must comply with the probability equation below:

ln[e K − Fi (e K − e−K )] ⎫ ϕi = ϕ¯ ± cos−1 ⎧ ⎬ ⎨ K ⎭ ⎩

T

∑ Psq sr (q−1) Dq = 1 s=1

(9)

(8) where K[−] is the Fisher constant; Fi[−] is a random number in the range of [0, 1]; and ϕ [deg ] is the mean value of the orientation angle. Fig. 1(a) illustrates a typical two-dimensional DFN model consisting of two fracture sets. The distribution parameters of the fracture networks are shown in Tables 1 and 2. The ratio of number densities between two fracture sets is 2 : 3.

where sr[−] is the length ratio between the side of a parent domain and its subdomain; Dq[−] is the multifractal dimension of centers at the dimension of order q[−]; and T is the total number of subdomains in the parent domain. This fractal distribution of location mimics the nature of fracture clustering. On the other hand, the Poisson process is the most commonly used assumption in DFN simulations, which distributes fracture centers randomly across the testing domain. Although less realistic than the multiplicative cascade process, the uniformly random distribution is easy to apply. The fracture orientation categorizes fractures into different sets. A 6

Engineering Geology 260 (2019) 105228

G. Ma, et al.

Fig. 6. Effects of fracture length distribution: (a) a = 2.0, p = 30; (b) a = 2.5, p = 30; (c) a = 2.0, p = 50; (d) a = 2.5, p = 50; (e) errors in permeability curves.

3. Equivalence of discrete fracture networks

retains most characteristics of the original domain. Hence, there exists an equivalent permeability factor that represents the proportion of transmissibility obtained before and after density reduction. Results in the original domain can be restored from the simulation of the equivalent model. In this paper, a correlation index S is defined to evaluate the significance of individual fractures in the overall fracture networks, which yields

3.1. Correlation index The fundamental concept of this equivalent method is that the contribution of individual fractures to the hydraulic properties in fracture networks is evaluable based on their geometric characteristics. The direct relations between geometric data (e.g., pore size and void ratio) and hydraulic conductivity have been widely accepted in wide range of sediments (Ren and Santamarina, 2018). Different from monoporous medium, hydraulic properties in fracture networks exhibits much more significant anisotropy and relates to more complex geometric parameters. This paper limits the scope of study to a parameter set including fracture lengths, apertures, and locations. Accordingly, the significance of each fracture is then utilized in fracture-density reduction. By neglecting less relevant fractures, the skeleton of networks

S=

N¯ b¯ γ P¯cβ

(10)

where N is the normalized intersection number that can be derived from the ratios of the intersection number in the current fracture to the average value of all fractures. The intersection number denotes the contribution of fractures to the domain connectivity, commonly proportional to its length. The normalized hydraulic aperture b 7

Engineering Geology 260 (2019) 105228

G. Ma, et al.

Fig. 7. Meshes generated in fracture-matrix domain with different extent of fracture reduction: (a) R0; (b) R50; (c) R70; (d) R90.

reduced domain. For examples, R0, R50, R70, and R90 indicate domains with 0%, 50%, 70%, and 90% of the total fracture neglected in sequence. As apertures are a constant in this case, the fracture correlation index depends on its locations and the number of intersections. When the first 50% of fractures are screened, the DFN configuration is scarcely changed since the most irrelevant half of fractures only take 36.68% of mass density and are poorly connected with minimal effects on the overall hydraulic properties. When the number density of fracture decreases by 70%, the mass density drops to 43.04%. Fractures in dense areas are primarily overlooked, and their influence on the local area and the entire domain can be substituted by the adjacent ones. When the density reduction reaches up to 90%, only 18.63% of mass density is retained. At this stage, the improvement in computational efficiency will sacrifice a massive amount of information on discontinuities, and the result accuracy will be doubtable.

Table 3 Results of fracture-matrix domain. Reduction level

E0 E10 E20 E30 E40 E50 E60 E70 E80 E90

Degree of freedom

Average aspect ratio

Error (εr)

[−]

[−]

[%]

20979 21015 20350 19624 18794 17780 16719 15656 14043 11811

0.74212 0.74275 0.75216 0.76451 0.77873 0.79715 0.81811 0.83919 0.87013 0.90692

0.00 0.00 0.09 0.32 0.75 1.17 1.45 2.38 5.32 7.92

additionally contributes to the fracture zone transmissibility. This contribution is then weighed by the factor γ. Based on the cubic law assumption, we use γ = 2 in two-dimensional fracture networks and γ = 3 for the three-dimensional space. In contrast, the importance of a fracture is reversely related to its density probability Pc since the fracture located at the denser area will have better substitutability. This density probability is also normalized into Pc and weighed by β. In this paper, we assume that β = 1 if the locations of fracture centers are generated from a multiplicative cascade process and β = 0 when the uniformly random distribution applies. Based on the correlation index, the fracture density can be reduced accordingly. Fig. 1 shows a typical fracture reduction process of the dense DFN model. Abbreviations represent the result of the density-

3.2. Equivalent permeability factor The tow-point flux approximation (TPFA) is utilized to generate the directional permeability curve that represents the transmissibility of a domain. Because an impermeable rock matrix is assumed, the flow field can be readily determined by a pipe-network method (Ren et al., 2015; Xu et al., 2018). For a static single-phase flow, mass conservation at intersections can be derived with the classic Darcy's law. Conventional methods represent the permeability anisotropy of fractured media through the ratio of the maximum and the minimum values among different directions (Barton and Quadros, 2015). And the permeability is described with a tensor, by which the directional permeability curve is an ellipse in two dimensions. However, the curve 8

Engineering Geology 260 (2019) 105228

G. Ma, et al.

Fig. 8. Contour of pressure in equivalent domains with different extent of fracture reduction: (a) E0; (b) E50; (c) E70; (d) E90.

hydraulic properties of the fractured domain. The enveloped area of this curve indicates the overall transmissibility, while the irregular perimeter shape represents the heterogeneity and anisotropy. During fracture density reduction, the domain transmissibility declines, while the loss of heterogeneity and anisotropy shows a lag effect. As shown in Fig. 3(a), the transmissibility in all directions is reduced when 50% of fractures are neglected. Yet, the permeability curve shows apparent shape similarity with results from the original domain despite the complex fracture configurations. This shape similarity represents the preservation of discontinuity. Accordingly, there exists an equivalent permeability factor ωr that can restore the original permeability curve from a density-reduced domain because of the shape similarity and equal proportion of transmissibility in all directions. The definition of equivalent permeability factor is expressed below (Ren et al., 2017).

seems more complex in a fractured system and cannot be described with a ratio or a tensor. In the following tests, we obtain the directional permeability curve by a rotational subdomain. As shown in Fig. 2, the subdomain simulated has a size of 5m × 5m, and the domain rotation angle is 15∘ at each step. Isolated fractures and fracture dead-ends are overlooked to obtain a backbone of fracture networks. With a constant pressure gradient imposed on the opposite boundaries, the average transmissibility of the domain can be derived in the testing directions. To focus on the pattern of permeability curve, the results are represented by the normalized values and expressed as follows:

k k¯rθ = rθ k 00

(11)

where krθ[L2] is the average permeability of domain along the testing direction. The subscripts, r and θ, represent the extent of density reduction and the testing direction, respectively. k00 represents the result when the constant pressure gradient is applied to the original domain along the x-axis. The unique pattern of the permeability curve illustrates the

Mr =

ωr = 9

1 2

∫0 M0 Mr

2π

2 k¯rθ dθ

(12)

(13)

Engineering Geology 260 (2019) 105228

G. Ma, et al.

εr =

1 2

2π ∫0 |(ωr k¯rθ )2 − k¯ 02θ | dθ

M0

(14)

This parameter is regarded as the error due to information loss and determines the validity of this equivalent method. As shown in Fig. 3, during density reduction, increasingly relevant fractures are neglected, and the accuracy of equivalent curves diminishes. When 50% of the fractures are neglected, the equivalent permeability curve perfectly agrees with the original one by an equivalent permeability factor ω50 down to 1.52, and the variance in similarity ε50 is as low as 0.060. When 70% of the fractures are overlooked, the permeability curve shrinks slightly more. In this case, the factor for restoration rises to 2.34, and the gap between the equivalent and original curves enlarges with an error up to 0.137. As the extent of density reduction approaches 90%, the resulting curve differentiates from the original one dramatically. To restore the permeability curve, the equivalent factor increases exponentially to 5.22, while the error induced soars to 0.280. This parameter also describes the result confidence of the simplified domain, and hence, it can be used in the limit state design (LSD) method. In practice, a lower bound exists in the number of the most critical and deterministic fractures that meets the requirements of precision and certainty. For example, if the limit of error based on design requirements is set at 0.100 in the demonstrated case, at least 32% of the most significant fractures should be included in the simulation. In addition to the extent of density reduction, the result accuracy also varies with different realizations based on the same statistics of fracture distribution. The certainty of the error range is a function of the reduction level and the distribution models. 4. Case study and discussion 4.1. Effects of fracture aperture distributions Fig. 9. Simulation domain setups of gas reservoir: (a) geometric configuration of reservoir; (b) boundary conditions and location of recovery well.

Two types of distribution patterns are considered in the following three cases to examine the fracture aperture effect. For the first, fracture apertures are related to their length since the in-situ stress tends to enhance the hydraulic aperture of the long fracture. The other two cases are designed to compare the effects of standard deviations when a lognormal distribution is applied to the fracture apertures. The standard deviations are assumed as 20% and 50% of the mean value, respectively. The detailed parameters of aperture distributions are listed as below:

Table 4 Normalized errors in estimation of gas production under different designs [%]. Reduction level

E0 E10 E20 E30 E40 E50 E60 E70 E80 E90 a

b

Wells 1 and 2b

Operations with single well Well 1a

Well 2a

Well 3a

0.00 0.10 0.00 0.08 0.22 0.33 0.16 0.70 1.70 1.14

0.00 0.10 0.20 0.73 1.14 2.62 3.02 4.98 8.60 15.31

0.00 0.12 0.08 0.25 0.53 0.90 0.87 2.17 4.44 5.13

0.00 0.27 0.21 0.34 0.12 0.91 1.67 1.43 0.30 8.36

KIC⋅=⋅0.6⋅MPa · m½, E = 5 × 104MPa, v = 0.25 for length-related apertures; μb = 5mm, σb = 1mm for log-normal distributed apertures with small deviation; μb = 5mm, σb = 2.5mm for log-normal distributed apertures with large deviation. The geometric configuration of fracture networks is identical with the previous case, which has a constant aperture distribution. Ten realizations are conducted to present stable results from log-normal distributions. As illustrated in Fig. 4, the variation in apertures will reshape the permeability curves, while the shape similarity still holds in each realization with the suggested equivalent method. Fig. 4(a) to (d) show the shape similarity in a typical realization of four different aperture distributions. Comparing the constant apertures in Fig. 4(e), if the fracture opening is related to its length, errors induced during density reduction are depressed since the significance of long fractures is enhanced and the skeleton of fracture networks becomes more relevant in respecting the whole domain. In contrast, if the fracture aperture follows the lognormal distribution, errors in the shape similarity of permeability curves increases at an accelerated rate given the deterioration of information loss. The average results from ten realizations are

Only a single borehole performs as the producing well. Wells 1 and 2 perform as the injection and production wells, respectively.

where M0 and Mr are the enveloped areas of the permeability curves in the original and density-reduced domains, respectively. Besides the transmissibility in the fluid flow, this equivalent permeability factor also applies to the diffusivity in the mass and thermal transport.

3.3. Similarity in permeability curves As the density reduction increases, the hydraulic characteristics all fade away. We define the variance in similarity εr to measure the deviation of hydraulic properties between the original and equivalent domains, which yields 10

Engineering Geology 260 (2019) 105228

G. Ma, et al.

connectivity are identified according to the length exponent. For a = 2.5, fractures smaller than the system size rule the domain's connectivity. Therefore, the reduction in fracture density will induce errors at the early stage, and the variance in simularity builds up gradually. As shown in Fig. 6(e), it reaches a limit of 10% when approximately 85% of fractures are neglected. The percolation coefficient has minor effects on the equivalent process in these two cases with a higher length exponent. In contrast, when a = 1.5, the connectivity is dominated by the largest fractures in the system. In these cases, the skeleton of fracture networks better represents the original system. Hence, disregarding the majority of fractures with short lengths has insignificant impact on the hydraulic properties of the domain. Even when the extent of reduction is 40%, the average errors in cases when p = 30 and p = 50 are only 0.025 and 0.039, respectively. As the reduction rate increases, the average errors in these two cases soar at an accelerated rate and exceed those in the cases with a higher length exponent.

approximately identical between two log-normal distributions at the same reduction rate. With a significant standard deviation in aperture distribution, the errors among ten realizations will be increasingly diverse. Thus, the results obtained from the equivalent method become less reliable. The distribution of correlation factors can explain the fracture aperture effects. Both constant apertures and log-normal distributions will reduce the significance of long fractures. This is especially noteworthy in log-normal distributions as the increasing randomness in apertures will narrow the difference of the correlation factor among individual fractures. These arguments will all weaken the applicability of this proposed method. 4.2. Effects of fracture location distributions In this study, the locations of fracture centers follow a fractal model and are generated by a multiplicative cascade process. The fractal dimension of centers determines the extent of clustering. When it equals to 2, the fractures are evenly distributed across the simulation domain. As the fractal dimension decreases, fractures tend to cluster, and the domain is more likely heterogeneous. Three cases with different fractal dimensions are designed and compared with the DFN models generated from a Poisson process. In each case, ten realizations are simulated to stabilize results. Although each realization has significant randomness, the similarity in the shape of permeability curves exists independently during the process of density reduction. As shown in Fig. 5(a) to (d), the typical configurations of fracture distribution are illustrated with different clustering extents. When the fracture networks are generated from a Poisson process, which is equivalent to the model with Dq = 2, the fracture centers tend to be uniformly distributed. In this case, the simulation domain performs as a homogeneous medium (i.e., β = 0). As the fracture density is reduced evenly across the domain, the errors induced remain stable below 0.082, as shown in Fig. 5(e). As the fractal dimension decreases, the mass density of fractures in different zones becomes more diverse. In these cases, the fracture correlation index shall include the effect of fracture number density (i.e., β = 1). Fractures in the denser area are neglected primarily since their contribution to the fluid flow is more substitutable. Consequently, errors are depressed at the early stage of fracture reduction, especially in the domains with a lower fractal dimension of location. As the fracture density is reduced further, the higher extent of fracture clustering will lead to more significant variance in the shape similarity of permeability curves. Considering εr = 0.1 as the upper limit of the error, this equivalent method is not applicable when the extents of fracture reduction reach 81.6%, 91.3%, and 98.0% for the cases with Dq = 1.75, Dq = 1.85, and Dq = 1.95, respectively. The standard deviation of results from ten realizations remains approximately identical at the same reduction rate despite the difference in distribution pattern.

5. Applications 5.1. The equivalent model for fracture-matrix domain With consideration for the rock matrix, reducing the fracture density with the equivalent model is even more significant. In this part, a fractured porous medium is generated containing 3, 340 discontinuities in a 5m × 5m domain. The upper and lower bounds of fracture lengths are set at 15m and 0.1m, respectively. The length exponent is 2.5, and the percolation coefficient is 20. In this case, fracture centers are uniformly distributed over the entire domain, and the rest of the parameters for DFN generation are identical to those in Tables 1 and 2. The permeability in the rock matrix is small enough to highlight the effects of fractures, and a constant pressure gradient is assigned between the lower left and upper right corners. Fig. 7 shows the resultant meshes at different extents of density reduction. The meshes are generated with a self-programmed software (Wang et al., 2017). Their quality is measured by the average aspect ratio (Hyman et al., 2014; Wang et al., 2017), and the total degree of freedom (i.e., the number of nodes) is used to represent the simulation domain complexity. As shown in Table 3, the total degree of freedom decreases with the number of fractures, while the quality of meshes increases since the fracture-density reduction will limit the number of ill-conditioned constraints. Even when 90% of fractures are neglected, the domain transmissibility in equivalence still has perfect agreement with the original one, and the error induced is as low as 7.92%. Moreover, the total degree of freedom is cut approximately in half, and the quality of meshed elements is elevated significantly. The effects of density reduction on the local pressure distribution are also compared in Fig. 8. The trend of the pressure head remains identical since the skeleton of fractures is retained in all levels of fracture reduction. However, the local details vanish gradually. When the reduction rate is 50%, the patterns of pressure contour are almost the same before and after the equivalent process. As 70% of fractures are neglected, only significant information is recorded in the equivalent domain. At a 90% reduction rate, most local details of the original domain are lost, and the skeleton of DFN cannot represent the original discontinuities.

4.3. Effects of fracture length distributions In a classic fractal model, the number density of fractures varies with the fractal dimension of length and the detectable minimum length. The mass density of fractures is also a function of the percolation coefficient. With constant upper and lower bounds of fracture lengths included in simulations, four cases are designed to compare the effects of length distributions. In these cases, the Poisson process is utilized to generate the fracture center locations, while other parameters of the DFN configurations remain identical to those in Tables 1 and 2. Ten realizations are simulated with each combination of length exponent and percolation coefficient to make the presentation of results quantitatively stable. Fig. 6(a) to (d) demonstrates the typical geometric configuration generated from the four designed cases. Larger percolation coefficients represent better connectivity in the system, while different regions of

5.2. Optimization of reservoir design The gas production in an NFR is simulated to validate a possible application of this equivalent model to a practical problem. The geometric configuration of the experiment model is illustrated in Fig. 9(a). A curving reservoir zone is embedded in the geologic formation, 1, 000m long and 800m deep. The upper and lower bounds of the reservoir zone can be represented as two sine functions. The parameters for DFN distributions are identical to those in Tables 1 and 2 excluding the maximum and minimum fracture length, which are set as 10m and 1, 11

Engineering Geology 260 (2019) 105228

G. Ma, et al.

350m, respectively. Several designs are conducted and compared based on different production wells to optimize the gas production in the reservoir. As shown in Fig. 9(b), a possible production well is drilled at x = 300m vertically with a depth of 500m (i.e., Well 1), the second design of a production well is located at x = 800m and y = 305m (i.e., Well 2), and a horizontal well (i.e., Well 3) is also considered, which has a 500m vertical segment at x = 500m, connected with a short curving part and a 150m horizontal segment. Accordingly, four designs are tested to estimate gas production. In the first three cases, only one single bore hole is activated in each design. A constant pressure head is applied on Boundaries 1 and 2, while Boundaries 3 and 4 are considered as nonflowing boundaries. The pressure head at the active well is regarded as 0. In the fourth case, all boundaries are impermeable, and Wells 1 and 2 are activated as the injection and production wells, respectively. Table 4 shows the normalized error induced during the process of density reduction in each individually designed case. In general, the errors are limited under 15.31% among all four cases, which makes this equivalent method applicable in the optimization of reservoir design. When the production well is located at the position of Well 1, the results from the density-reduced domain are in accordance with the accurate result. The maximum deviation in productivity is as low as 1.14%, even when 90% of fractures are neglected. In contrast, when the production well is located at the position of Well 2, errors at the same extent of fracture reduction become dramatic, which limits the maximum fracture density reduction rate. When the horizontal well is used, the gap between equivalent and original production is still manageable. Similarly, if Wells 1 and 2 perform as the injection and production wells, respectively, the maximum error in all tested equivalent domains is limited under 10%. Therefore, the estimation obtained from densityreduced DFN models with a selected fracture skeleton has excellent agreement with the accurate result, which can be adapted in the optimization of reservoir designs.

•

reduces computational complexity but also preserves overall transmissibility and local details of distorted pressure field with negligible errors. We also apply the density reduction scheme to fast and robust production estimations in the optimization of gas reservoir designs. It is noted that the proposed approach is preferred in reservoirs with fracture clustering, in which the fracture networks cannot be treated as an equivalent continuum.

We should emphasize that the accuracy of our equivalent model is highly dependent on the accessibility of site information in engineering applications. For the time being, characterizing distribution models and regenerating discrete fracture networks from site data are frequently difficult if not impossible. Another challenge for the equivalent model is the certainty analysis of simulation results. Results uncertainty is not only induced among various iterations of fracture networks regenerations but also deteriorated by fracture reduction processes. The effects of result uncertainty on the computational complexity is not included in this study. Furthermore, in terms of mechanical behavior, the significance of individual fractures is crucial but still inconclusive. Hence, future works would extend the correlation index to evaluate effects of fracture on rock strength and deformability. In addition, error induced by fracture density reduction remains inevitable. Implementations that further improve the model accuracy and result confidence need to be discussed. Acknowledgment This research was supported by the Research Training Program (RTP) scheme at the University of Western Australia (UWA). The authors wish to thank the National Natural Science Foundation of China (NSFC) for their financial support (No. 51778029 & No. 51627812). We would like to thank the reviewers for their constructive suggestions and discussion.

6. Conclusion References This paper proposes an improved equivalent DFN model to generate a skeleton of fracture networks that can sufficiently represent the original domain. The fundamental concept of this study is to quantitatively analyze the contribution of individual fractures on the hydraulic properties of a fractured medium. We define the correlation factor to determine the significance of fractures, from which the density-reduced fracture networks are derived.The fracture skeleton retains most characteristics of the original domain and can restore its transmissibility by a simple parameter of equivalent permeability factor. The equivalent model improves the computational efficiency by reducing the number of fractures and depressing the total degree of freedom in the simulation domain. Furthermore, the sparser fracture density will ease the computational burden in the discretization of a fracture-matrix system and elevate the quality of meshed elements through the reduction in illconditioned intersections. Since the fracture skeleton is preserved, the connectivity and transmissibility of the original domain can be conveniently restored with a manageable level of error and uncertainty. In this paper, the following works have been done:

Baghbanan, A., Jing, L., 2007. Hydraulic properties of fractured rock masses with correlated fracture length and aperture. Int. J. Rock Mech. Min. Sci. 44, 704–719. https://doi.org/10.1016/J.IJRMMS.2006.11.001. Barton, N., Quadros, E., 2015. Anisotropy is everywhere, to see, to measure, and to model. Rock Mech. Rock. Eng. 48, 1323–1339. https://doi.org/10.1007/s00603-0140632-7. Bear, J., 1988. Dynamics of Fluids in Porous Media. (Dover, New York). Berkowitz, B., 1995. Analysis of fracture network connectivity using percolation theory. Math. Geol. 27, 467–483. https://doi.org/10.1007/BF02084422. Bidgoli, M.N., Zhao, Z., Jing, L., 2013. Numerical evaluation of strength and deformability of fractured rocks. J. Rock Mech. Geotech. Eng. 5, 419–430. https://doi.org/ 10.1016/J.JRMGE.2013.09.002. Bour, O., Davy, P., 1997. Connectivity of random fault networks following a power law fault length distribution. Water Resour. Res. 33, 1567–1583. https://doi.org/10. 1029/96WR00433. Chen, J., Lan, H., Macciotta, R., Wu, Y., Li, Q., Zhao, X., 2018. Anisotropy rather than transverse isotropy in Longmaxi shale and the potential role of tectonic stress. Eng. Geol. 247, 38–47. https://doi.org/10.1016/J.ENGGEO.2018.10.018. Chen, Y., Ma, G., Li, T., Wang, Y., Ren, F., 2018a. Simulation of wormhole propagation in fractured carbonate rocks with unified pipe-network method. Comput. Geotech. 98, 58–68. https://doi.org/10.1016/j.compgeo.2017.11.009. Chen, Y., Ma, G., Wang, H., Li, T., 2018b. Evaluation of geothermal development in fractured hot dry rock based on three dimensional unified pipe-network method. Appl. Therm. Eng. 136, 219–228. https://doi.org/10.1016/J.APPLTHERMALENG. 2018.03.008. Chen, Y., Ma, G., Jin, Y., Wang, H., Wang, Y., 2019a. Productivity evaluation of unconventional reservoir development with three-dimensional fracture networks. Fuel 244, 304–313. https://doi.org/10.1016/J.FUEL.2019.01.188. Chen, Y., Ma, G., Wang, H., Li, T., Wang, Y., 2019b. Application of carbon dioxide as working fluid in geothermal development considering a complex fractured system. Energy Convers. Manag. 180, 1055–1067. https://doi.org/10.1016/J.ENCONMAN. 2018.11.046. Chesnaux, R., Allen, D.M., Jenni, S., 2009. Regional fracture network permeability using outcrop scale measurements. Eng. Geol. 108, 259–271. https://doi.org/10.1016/J. ENGGEO.2009.06.024. Darcel, C., Bour, O., Davy, P., de Dreuzy, J.R., 2003. Connectivity properties of twodimensional fracture networks with stochastic fractal correlation. Water Resour. Res.

• Stochastic DFN models are generated in two dimensions based on a • • •

set of fracture-distribution parameters of length, location, orientation, and aperture. A case is investigated to assess the similarity of permeability curves generated from equivalent models with different fracture density, and thereafter, the result accuracy is validated. The effects of fracture configurations on the applicability of the equivalent model are also studied by adjusting the parameters of fracture distributions in DFN generation. This model is adapted in the preprocessing of mesh generation in a fracture-matrix domain. The fracture skeleton not only significantly 12

Engineering Geology 260 (2019) 105228

G. Ma, et al.

Maryška, J., Severýn, O., Vohralík, M., 2005. Numerical simulation of fracture flow with a mixed-hybrid FEM stochastic discrete fracture network model. Comput. Geosci. 8, 217–234. https://doi.org/10.1007/s10596-005-0152-3. Meakin, P., 1991. Invasion percolation on substrates with correlated disorder. Physica 173, 305–324. https://doi.org/10.1016/0378-4371(91)90366-K. Olson, J.E., 2003. Sublinear scaling of fracture aperture versus length: an exception or the rule? J. Geophys. Res. Solid Earth 108. https://doi.org/10.1029/2001JB000419. Pan, J.B., Lee, C.C., Lee, C.H., Yeh, H.F., Lin, H.I., 2010. Application of fracture network model with crack permeability tensor on flow and transport in fractured rock. Eng. Geol. 116, 166–177. https://doi.org/10.1016/J.ENGGEO.2010.08.007. Parashar, R., Reeves, D.M., 2012. On iterative techniques for computing flow in large two-dimensional discrete fracture networks. J. Comput. Appl. Math. 236, 4712–4724. https://doi.org/10.1016/J.CAM.2012.02.038. Priest, S.D., 1993. Discontinuity Analysis for Rock Engineering. Springer, Netherlands. Reeves, D.M., Parashar, R., Pohll, G., Carroll, R., Badger, T., Willoughby, K., 2013. The use of discrete fracture network simulations in the design of horizontal hillslope drainage networks in fractured rock. Eng. Geol. 163, 132–143. https://doi.org/10. 1016/J.ENGGEO.2013.05.013. Ren, X.W., Santamarina, J.C., 2018. The hydraulic conductivity of sediments: a pore size perspective. Eng. Geol. 233, 48–54. https://doi.org/10.1016/J.ENGGEO.2017.11. 022. Ren, F., Ma, G., Fu, G., Zhang, K., 2015. Investigation of the permeability anisotropy of 2D fractured rock masses. Eng. Geol. 196, 171–182. https://doi.org/10.1016/J. ENGGEO.2015.07.021. Ren, F., Ma, G., Fan, L., Wang, Y., Zhu, H., 2017. Equivalent discrete fracture networks for modelling fluid flow in highly fractured rock mass. Eng. Geol. 229, 21–30. https:// doi.org/10.1016/J.ENGGEO.2017.09.013. Renshaw, C.E., 1996. Influence of subcritical fracture growth on the connectivity of fracture networks. Water Resour. Res. 32, 1519–1530. https://doi.org/10.1029/ 96WR00711. Renshaw, C.E., Park, J.C., 1997. Effect of mechanical interactions on the scaling of fracture length and aperture. Nature 386, 482–484. https://doi.org/10.1038/ 386482a0. Schertzer, D., Lovejoy, S., 1987. Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes. J. Geophys. Res. 92, 9693. https://doi. org/10.1029/JD092iD08p09693. Vu-Bac, N., Lahmer, T., Zhuang, X., Nguyen-Thoi, T., Rabczuk, T., 2016. A software framework for probabilistic sensitivity analysis for computationally expensive models. Adv. Eng. Softw. 100, 19–31. https://doi.org/10.1016/J.ADVENGSOFT. 2016.06.005. Wang, Y., Ma, G., Ren, F., Li, T., 2017. A constrained Delaunay discretization method for adaptively meshing highly discontinuous geological media. Comput. Geosci. 109, 134–148. https://doi.org/10.1016/J.CAGEO.2017.07.010. Xu, C., Dowd, P.A., Mardia, K.V., Fowell, R.J., 2007. A connectivity index for discrete fracture networks. Math. Geol. 38, 611–634. https://doi.org/10.1007/s11004-0069029-9. Xu, C., Fidelibus, C., Dowd, P., Wang, Z., Tian, Z., 2018. An iterative procedure for the simulation of the steady-state fluid flow in rock fracture networks. Eng. Geol. 242, 160–168. https://doi.org/10.1016/J.ENGGEO.2018.06.005. Yang, T., Liu, H.Y., Tang, C.A., 2017. Scale effect in macroscopic permeability of jointed rock mass using a coupled stress–damage–flow method. Eng. Geol. 228, 121–136. https://doi.org/10.1016/J.ENGGEO.2017.07.009. Yin, Q., Ma, G., Jing, H., Wang, H., Su, H., Wang, Y., Liu, R., 2017. Hydraulic properties of 3D rough-walled fractures during shearing: an experimental study. J. Hydrol. 555, 169–184. https://doi.org/10.1016/J.JHYDROL.2017.10.019. Yu, B., Zou, M., Feng, Y., 2005. Permeability of fractal porous media by Monte Carlo simulations. Int. J. Heat Mass Transf. 48, 2787–2794. https://doi.org/10.1016/J. IJHEATMASSTRANSFER.2005.02.008. Zhang, C., Chen, Q., Qin, X., Hong, B., Meng, W., Zhang, Q., 2017. In-situ stress and fracture characterization of a candidate repository for spent nuclear fuel in Gansu, northwestern China. Eng. Geol. 231, 218–229. https://doi.org/10.1016/J.ENGGEO. 2017.10.007. Zhao, Y., Feng, Z., Liang, W., Yang, D., Hu, Y., Kang, T., 2009. Investigation of fractal distribution law for the trace number of random and grouped fractures in a geological mass. Eng. Geol. 109, 224–229. https://doi.org/10.1016/J.ENGGEO.2009.08.002. Zhou, S., Zhuang, X., Rabczuk, T., 2018. A phase-field modeling approach of fracture propagation in poroelastic media. Eng. Geol. 240, 189–203. https://doi.org/10. 1016/J.ENGGEO.2018.04.008.

39, 1272. https://doi.org/10.1029/2002WR001628. de Dreuzy, J.R., Davy, P., Bour, O., 2001. Hydraulic properties of two-dimensional random fracture networks following a power law length distribution: 1. Effective connectivity. Water Resour. Res. 37, 2065–2078. https://doi.org/10.1029/ 2001WR900011. de Dreuzy, J.R., Davy, P., Erhel, J., de Brémond d’Ars, J., 2004. Anomalous diffusion exponents in continuous two-dimensional multifractal media. Phys. Rev. E 70, 016306. https://doi.org/10.1103/PhysRevE.70.016306. Earon, R., Olofsson, B., 2018. Hydraulic heterogeneity and its impact on kinematic porosity in Swedish coastal terrains. Eng. Geol. 245, 61–71. https://doi.org/10.1016/J. ENGGEO.2018.08.008. Hamdia, K.M., Silani, M., Zhuang, X., He, P., Rabczuk, T., 2017. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. Int. J. Fract. 206, 215–227. https://doi.org/10.1007/s10704-0170210-6. Hamdia, K.M., Ghasemi, H., Zhuang, X., Alajlan, N., Rabczuk, T., 2018. Sensitivity and uncertainty analysis for flexoelectric nanostructures. Comput. Methods Appl. Mech. Eng. 337, 95–109. https://doi.org/10.1016/J.CMA.2018.03.016. Hamm, S.Y., Kim, M., Cheong, J.Y., Kim, J.Y., Son, M., Kim, T.W., 2007. Relationship between hydraulic conductivity and fracture properties estimated from packer tests and borehole data in a fractured granite. Eng. Geol. 92, 73–87. https://doi.org/10. 1016/J.ENGGEO.2007.03.010. Huang, N., Jiang, Y., Liu, R., Li, B., 2017. Estimation of permeability of 3-D discrete fracture networks: an alternative possibility based on trace map analysis. Eng. Geol. 226, 12–19. https://doi.org/10.1016/J.ENGGEO.2017.05.005. Hyman, J.D., Gable, C.W., Painter, S.L., Makedonska, N., 2014. Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy. SIAM J. Sci. Comput. 36, A1871–A1894. https://doi.org/10.1137/130942541. Ko, N.Y., Ji, S.H., Koh, Y.K., Choi, J.W., 2015. Evaluation of two conceptual approaches for groundwater flow simulation for a rock domain at the block-scale for the Olkiluoto site, Finland. Eng. Geol. 193, 297–304. https://doi.org/10.1016/J. ENGGEO.2015.05.003. Larsson, M., Odén, M., Niemi, A., Neretnieks, I., Tsang, C.F., 2013. A new approach to account for fracture aperture variability when modeling solute transport in fracture networks. Water Resour. Res. 49, 2241–2252. https://doi.org/10.1002/wrcr.20130. Le Borgne, T., Bour, O., Paillet, F.L., Caudal, J.P., 2006. Assessment of preferential flow path connectivity and hydraulic properties at single-borehole and cross-borehole scales in a fractured aquifer. J. Hydrol. 328, 347–359. https://doi.org/10.1016/J. JHYDROL.2005.12.029. Lee, C.H., Deng, B.W., Chang, J.L., 1995. A continuum approach for estimating permeability in naturally fractured rocks. Eng. Geol. 39, 71–85. https://doi.org/10.1016/ 0013-7952(94)00064-9. Li, B., Liu, R., Jiang, Y., 2016. Influences of hydraulic gradient, surface roughness, intersecting angle, and scale effect on nonlinear flow behavior at single fracture intersections. J. Hydrol. 538, 440–453. https://doi.org/10.1016/J.JHYDROL.2016.04. 053. Liu, R., Jiang, Y., Li, B., Wang, X., 2015. A fractal model for characterizing fluid flow in fractured rock masses based on randomly distributed rock fracture networks. Comput. Geotech. 65, 45–55. https://doi.org/10.1016/J.COMPGEO.2014.11.004. Liu, R., Li, B., Jiang, Y., 2016. A fractal model based on a new governing equation of fluid flow in fractures for characterizing hydraulic properties of rock fracture networks. Comput. Geotech. 75, 57–68. https://doi.org/10.1016/J.COMPGEO.2016.01.025. Liu, R., Zhu, T., Jiang, Y., Li, B., Yu, L., Du, Y., Wang, Y., 2018. A predictive model correlating permeability to two-dimensional fracture network parameters. Bull. Eng. Geol. Environ. 1–17. https://doi.org/10.1007/s10064-018-1231-8. Lu, Y., Liu, S., Weng, L., Wang, L., Li, Z., Xu, L., 2016. Fractal analysis of cracking in a clayey soil under freeze–thaw cycles. Eng. Geol. 208, 93–99. https://doi.org/10. 1016/J.ENGGEO.2016.04.023. Ma, G.W., Wang, H.D., Fan, L.F., Wang, B., 2017. Simulation of two-phase flow in horizontal fracture networks with numerical manifold method. Adv. Water Resour. 108, 293–309. https://doi.org/10.1016/J.ADVWATRES.2017.08.013. Ma, G., Li, T., Wang, Y., Chen, Y., 2019a. Numerical simulations of nuclide migration in highly fractured rock masses by the unified pipe-network method. Comput. Geotech. 111, 261–276. https://doi.org/10.1016/J.COMPGEO.2019.03.024. Ma, G., Wang, Y., Li, T., Chen, Y., 2019b. A mesh mapping method for simulating stressdependent permeability of three-dimensional discrete fracture networks in rocks. Comput. Geotech. 108, 95–106. https://doi.org/10.1016/J.COMPGEO.2018.12.016.

13