## Abstract

A methodology for simulating groundwater flow in three-dimensional (3D) stochastic fracture rocks based on a commonly used finite-difference method is presented in this paper. Different realizations of fracture networks are generated by the fracture continuum method (FCM), in which appropriate 3D cuboids are used to describe the geometry of fractures. Then, the effects of different parameter distributions on the fracture networks indicated that the length, orientation, and density of fractures all play significant roles in the connectivity of fractures in this methodology. Greater length and density and wider orientation range of fractures lead to greater connectivity. The proper contrast in hydraulic conductivities between the fractures and matrix is found to be approximately 10^{5} due to the contribution of fluid flow in the matrix which can be ignored. It is shown that the fracture density plays a key role in stabilizing the equivalent hydraulic conductivity ($Ke$) of the fracture networks. Furthermore, the greater length and closer orientation of fractures to the general flow direction, the larger $Ke$ of the generated fracture networks possess. The findings of this study can help for a better understanding of the mechanism of FCM and the influence of geometry characteristics on the hydraulic conductivity of FCM models.

## 1. Introduction

Fractured geologic formations are the repositories and flow channels of subsurface fluids (e.g., groundwater and petroleum). With the rapid development of industry, such formations are increasingly used as repositories to the dispose hazardous wastes, especially for the nuclear waste disposal sites. Therefore, the studies focusing on fractured geologic formations are frequently implemented to solve the subsurface flow problems, such as petroleum production, shale gas development, pit water management, and hazardous waste treatment [1–5]. Fractures are the main flow pathways for fluid flow and solute transport in fractured geologic formations. Thus, characterizing the fracture network and then modeling groundwater flow and transport in fractured media are definitely the vital issues associated with fractured geologic formations [6–9].

Based on the actual demand mentioned above, the methodology for simulating the flow and transport in fractured rocks, composed of the discrete fracture network (DFN) [10] and the stochastic continuum (SC) [8], has been developed. The equivalent discrete fracture approach was firstly designed on the hypothesis that fractures as the homogeneous porous medium by Long et al. [11]. Then, DFN approach was introduced by Dershowitz et al. [10] to model groundwater flow and solute transport in the complex fracture network, which is surrounded by impermeable rock matrix. The discrete fracture network was mainly applied to consider fracture networks as interconnected line segments in the two dimensions [12–17] and elliptic or polygonal discs in the three dimensions [18–23]. With the ability of accurate characterization and conceptualization, DFN requires the details of fracture network property such as the lengths and transmissivity, which is computationally expensive to achieve. This computational constraint restricts the application of DFN leading to smaller simulation scales (<100 m) [24]. Although the calculation ability of computer has been highly improved today, the challenges with a large number of input parameters, e.g., the location, orientation, shape, size, and transmissivity of each fracture in fractured media, further limit the applicability of the DFN at regional scale.

As an alternative to the DFN, the SC addresses the fractured aquifer system with an equivalent continuum. The approach is used to model groundwater flow and transport in the porous media based on the concept of the representative elementary volume (REV) [7]. SC assumes that, over some REVs, the fractured network can be characterized as the equivalent porous medium. Consequently, the mathematical model for groundwater flow in the equivalent porous media utilized the governing differential equation derived from Darcy’s law. The hydraulic conductivity can be considered as a regionalized variable following the random spatial distribution of block/element equivalent hydraulic conductivity that represents the spatially averaged fracture properties. One can use Monte Carlo method to generate stochastic conditional or unconditional realizations of hydraulic conductivity based on the analysis of the fracture systems. The SC approach, affording flexibility and ability to provide quantified ranges of behavior [25], can avoid the challenge for the detailed information on fracture geometry and assumptions associated with the nature of groundwater flow within individual fractures. However, in the practical applications, the SC cannot achieve accurate simulation results for the effect of individual fractures on fluid flow and solute transport.

The integrated framework, the fracture continuum method (FCM), combines the advantages of DFN and SC approach. In the FCM model, a fractured aquifer is discretely incorporated into a continuum model of groundwater flow by assigning high hydraulic conductivity to model grids occupied by the fracture network. The geometry of fractures is directly converted into the finite-difference grid to simulate groundwater flow [26, 27]. The FCM model incorporates the advantages of the DFN model, as well as an effective continuum representation of fracture conductivity with high processing speed [28]. Similar to the DFN method, if the location of a fracture is known, the fracture can be included in the FCM model. Otherwise, fractures are generated stochastically based on the probability distributions of fracture orientation and size derived from field observation and measurement. The flow in fractured reservoirs is mainly affected by the hydraulic gradient and the permeability of the fracture network. The hydraulic gradient is associated with the boundary conditions, while the permeability of the fracture network is dominated by the connectivity of the fracture network and the hydraulic conductivity of the fracture and matrix. Then, the overall permeabilities of fractured aquifers with different lithologic matrix of the bedrocks are different from each other. For example, illustrated in Figure 1 are the outcrops of two fractured aquifers with well-developed fracture networks in Anhui Province of China, where there are similar fracture apertures of approximately 1-10 cm, and the measured fracture dip angles are about 80 degrees in Figure 1(a) and 75 degrees in Figure 1(b), respectively. However, these two fractured aquifers possess different lithologic matrix: shown in Figure 1(a) is the fractured limestone beds in Tongling area, and shown in Figure 1(b) is the fractured argillaceous sandstone beds in Xuancheng area. However, due to the weaker permeability of argillaceous sandstone matrix compared to the limestone matrix, the difference between these two aquifers are obvious in water abundance, indicating that the single well capacity is only 50 m^{3}/d in the argillaceous sandstones, and is up to 500 m^{3}/d in the fractured limestones. In the DFN method, the matrix is generally considered to be impermeable, while the SC method does not make a special distinction between the fracture and the matrix. The FCM method can overcome this shortcoming. The generated geometry is used for mapping fracture network onto continuum grids, and differentiated permeability is used to represent fracture and matrix systems. In general, on the basis of retaining the geometric characteristics of the single fracture, the FCM method reduces the computational cost and preserves the connectivity between the fracture and matrix.

At present, most of the research on FCM methods focuses on 2D conditions. In these researches, fractures are converted to linear fissures in continuum models [29]. Successful mapping of discrete fractures onto regular grids provides an efficient and relatively simple way for simulating flow in fractured geologic media and allows the use of the traditional flow simulators for 2D solution in porous media [30]. In 3D conditions, the FCM model has been presented in heat transport in fractured systems [31, 32], where the permeability tensor has been applied to define permeability of each grid cell in the model [33]. However, the extended 3D models are still ambiguous for the solutions related to groundwater flow and transport, limiting their applications in real 3D fractured geological formations. Systematic studies on fracture mapping methods and hydraulic properties of fracture networks of 3D FCM models were rarely presented in field-scale real-world applications. To address this issue, this study attempts to construct a 3D FCM model that can efficiently translate discrete fractures generated from stochastic process onto 3D finite different grids. To evaluate the applicability of the 3D FCM model, the quality of generated stochastic fracture networks was examined to check the connectivity of fracture network in different cases of the model and investigate the basic characteristics of groundwater flow in the fracture-matrix field through test cases of stochastic 3D fracture networks. Figure 2 indicates the general sketch of the problem under this study. The stochastic fracture network generation algorithm and continuum mapping approaches are combined for transferring discrete fractures into FCM models. After the models are established, this study will focus on the flow simulation. All the algorithm used in this study was compiled by python language.

The remainder of this paper is organized as follows. First, the methodology is presented, in which the Geostatistical and Monte Carlo methods are combined to generate 3D fracture networks with prescribed geometry characteristics of fractures, including orientation, length, width, height, and density of fractures. Then, the generated fracture network is mapped onto 3D continuum grids according to the corresponding geometry characteristics. To eliminate unconnected grids, the percolation algorithm is proposed, and the results of different cases of fracture networks are presented. Geometry characteristics of fractures and the density of fractures in domain will be analyzed to inspect the quality and connectivity of the mapping method. On the basis of different cases of fracture networks generated, the contrast of hydraulic conductivity between different fractures and matrixes is also considered in different ratios. Furthermore, the comparison of effective conductivity among cases is discussed, and different conductivity ratios of fracture to rock matrix are investigated. When the effective ratio is determined, the flow solutions of different cases will be discussed to analyze the influence of fracture geometry characteristics on the equivalent hydraulic conductivity of the fracture network. Finally, the conclusions are summarized.

## 2. Methodology

### 2.1. Generation of Fracture Networks

The algorithm of generating stochastic individual fracture networks is learned from the FracSim3D code [34]. Distribution of fracture location points represented by the center points of fractures can be generated according to the homogenous Poisson process [35]. In the process, set $A$ was assumed as the family of subsets of region ($R$) of interest, while $NA$ is defined as the number of points in set $A$. Under the condition of $NA$ followed Poisson distribution ($\lambda vA$) for each $A$, in which $\lambda $ is the characteristic parameter called intensity and $vA$ is the volumetric measure of $A$, $NA1,NA2,\cdots ,NAn$ are independent for every finite collection $A1,A2,\cdots ,An$ of disjoint subsets, and $NA$ is a homogeneous Poisson point process with intensity $\lambda >0$. The homogenous Poisson point process $NA$ is characterized by two fundamental properties that have already appeared as asymptotic ones: (1) the random number of points $NA$ in a bounded Borel set $A$ follows a Poisson distribution of mean $\lambda vA$ for constant $\lambda $ (intensity), where $vA$ means the volumetric measure of $A$; (2) for every finite collection $A1,A2,\cdots ,An$ of disjoint subsets, $NA1,NA2,\cdots ,NAn$ are independent.

The location of fracture location point combined with fracture geometry properties determines the spatial morphological feature of a certain fracture. Once the configuration of a set of fracture location points is generated, the random parameters with respective probability density functions that define the geometric properties of fractures are assigned to each point. Then, Monte Carlo simulations are used to produce the realizations of fracture network according to the distributions of defined parameters of fractures. Instead of ellipsoids [36], the cuboids are applied to define the fracture shapes generated in this study. Three sets of size parameters of fractures including length, height, and width are assigned to each point. For a certain 3D fracture, three angle parameters including dip angle ($\alpha $), strike angle ($\beta $), and rotation angle ($\gamma $) are needed for precisely mapping the orientation of the fracture in the Cartesian coordinate system (Figure 3).

### 2.2. Continuum Mapping Approach for Fracture Networks

Spatial length, height, and aperture are the geometry parameters of fractures with the shape of cuboids. Aperture is the opening width of an individual fracture. With the application of FCM model, fractures are idealized as blocks with dimensions in terms of length, width, and height [25, 26], in which the fracture aperture is represented by the width and the fracture width is expressed by the height. Fractures generated by stochastic method are mapped onto 3D continuum grids, and each fracture block occupies a certain amount of specific grids. As shown in Figure 4, the fracture block is deterministically mapped onto discretized finite-difference grids according to the grids volume occupied. During the mapping process, fracture-occupied grids (dark-filled) are assigned as hydraulic properties which represent the fracture hydraulic properties, while nonfracture-occupied grids (no-color-filled) are assigned as hydraulic properties of the rock matrix.

In a 2D FCM model, flow is restricted in a single pathway through fractures. This leads to an easy adjustment of the conductivity of each grid so that the proper discharge can be obtained in each fracture. In 3D situations, by contrast, flow in each fracture is much more complex. The concept of hydraulic conductivity tensor for each grid cell in the model domain needs to be used to define the conductivity of each fracture [36]. This study is focused on the performance of the connectivity caused by geometry of fractures and the equivalent hydraulic characteristics of the fracture network domain. Thus, the adjustment of the conductivity of each grid is not discussed in this study. Based on the method of assigning conductivity parameters in heterogeneous aquifer systems containing narrow connected pathways [37, 38], the matrix conductivity ($Km$) is defined as the hydraulic conductivity of rock matrix in each model, and the fracture conductivity ($Kf$) is defined as the hydraulic conductivity of fracture. Also, the ratio ($Kf/Km$) of fracture conductivity to matrix conductivity is specified to ensure the majority of flow will move through fracture grids in flow simulation.

The connectivity of fracture network is controlled by the degree of fracture interconnectivity. Isolated fracture clusters that do not contribute to connectivity and flow in the model may exist in the generated fracture network. Thus, it is necessary to eliminate the isolated fracture segments to avoid calculating deviation and increasing the computing efficiency. Compared with the DFN model in that only connected fractures contribute to flow, the FCM model is a continuum one, in which isolated fracture clusters and dead ends still contribute to total flow. Then, mapping fracture network onto continuum grids can also lead to an enhancement of connectivity in the model [30]. Therefore, to resolve this enhancement issue, a percolation algorithm is used for mapping the generated stochastic fracture networks onto finite-difference grids and removing the isolated fracture segments. The backbone of the fracture network is then presented in the discretized finite-difference grids.

### 2.3. Simulation of Groundwater Flow in Facture Networks

## 3. Application to Stochastic 3D Fracture Networks

### 3.1. Description of the Fracture Networks

In this study, a synthetic cuboid with the size of $200\u2009m\xd7100\u2009m\xd750\u2009m$ is used as model domain, and the refined grids of $1\u2009m\xd71\u2009m\xd71\u2009m$ are uniformly discretized for the domain of interest. Within the domain, two sets of fractures with the same properties except for orientations are generated. The orientation of fractures in each set is uniformly and randomly distributed in defined intervals. The number of fracture centers per unit area and the location of fracture centers follow the Poisson point process. Fracture height is set as a fixed value of 4 m, and fracture width is defined as the unit size of the grid. Fracture length in this study follows uniform and lognormal distribution with different mean of lengths ($lmean$). The details of distribution and geometric parameters used in this study are listed in Table 1. A total of five different cases is considered in this study. For each of the five different cases, 200 Monte Carlo realizations of fracture networks are generated. Cases 1 and 2 explore the effect of fracture length, while case 3 is of the same parameters as cases 1 and 2 except for the larger fracture density. The effect of variability in orientation is investigated for case 4. Case 4 follows the uniform distribution in (0, 360), while for other four cases (cases 1–3 and 5), the ranges of orientation are limited from 60 to 80 and from 100 to 120. Case 5 uses the lognormal distribution for fracture length, with the mean of $\mu =3.8$ and the standard deviation of $\sigma =0.4$, respectively.

Once the parameters of fracture networks have been determined, fractures can be mapped into the continuum models and truncated at the model boundaries, with the isolated ones eliminated from the original fracture networks. Fracture density is always expressed as linear frequency, trace length density, or fracture area density. For the continuum grid composed model, an index of the fracture volume occupied ratio was used, i.e., the sum of the fracture occupied grid volume divided by the total domain volume, to represent the fracture density ($\delta $) in the domain.

Determined by dividing the standard deviation of $Ke$, the uncertainty for the typical behavior of $Ke$ can be evaluated. In the same case, the relationship between the fracture volume occupied ratio and the equivalent hydraulic conductivity of the whole fracture network can be quantified by correlation analysis. Then, different cases will be compared to evaluate the relationship between conductivity and properties of the fracture network.

### 3.2. Results and Discussion

The spatial structures of fracture networks for one realization in different cases generated by 3D FCM are shown in Figures 6(a)–6(j). The fracture-occupied grids in the domain are explicitly displayed in dark, which are assigned hydraulic properties of fracture medium, while the rest of nonfracture-occupied grids are treated as matrix in white. Figures 6(a), 6(c), 6(e), 6(g), and 6(i) depict the original structures obtained from the parameters mentioned in Table 1, while Figures 6(b), 6(d), 6(f), 6(h), and 6(j) exhibit the backbone structures of fracture networks without isolated fracture segments for cases 1-5, respectively. Fractures in all 5 cases can be accurately translated onto the continuum grids, and then, the backbones and their connectivity between fractures will be used in the following flow simulations.

Due to the truncation and elimination of isolated fractures, the density of fracture networks becomes smaller than the original construction. For comparison, Figure 7 shows the ratios of the backbones to the original structures of multiple realizations of the generated fracture networks for the different cases 1-5. The influence of fracture geometry on the fracture network connectivity can be examined in Figure 7, in which a larger ratio means that more fractures are retained for later flow simulation, indicating better connectivity in that case. Cases 1 and 2 indicate that the fracture length is an important parameter to define connectivity of fracture networks with similar densities of fractures. A comparison of cases 1 and 2 also indicates that the portion of fracture segments truncated from the network decreases apparently as the fracture length increases. Further, the fracture volume occupied ratio for case 3 is much larger than that in other cases, and the portion of fractures left in backbones is similarly bigger than in other cases. It is due to that the higher fracture volume occupied ratio has a larger opportunity for the fractures to intersect each other in the fracture network. A greater ratio of the backbone and original structure for case 4 in comparison with case 1 as in Figure 7 indicates that orientation exerts influence on the connectivity of fracture networks. Moreover, a wider range of fracture orientation distribution always leads to more fracture intersections. Case 5 is similar to case 1 except that the fracture length for case 5 follows the lognormal distribution with a standard deviation of $\sigma $ equal to 0.4. However, there is only slight difference between lognormal and uniform distributions with the same average fracture length for cases 1 and 5.

As the ratio of $Kf$ to $Km\u2009Kf/Km$ is proposed in flow simulation, different combinations of $Kf$ and $Km$ can be used in simulation to check if the flow mainly occurs in fractures of fractured rocks. In this study, it is assumed that $Qtotal$ represents the total flow discharge (volume per unit time) moving from upstream to downstream boundaries. A total of 200 realizations of the generated stochastic fracture networks are then used for the analysis of case 1. For this purpose, a fixed value of $1.0\xd710\u22122\u2009m/s$ is assigned for $Kf$, with the value of $Kf/Km$ ranging from 10^{2} to 10^{6}. As shown in Figure 8(a), there is no evident difference in total flow between the generated fracture networks with the value of $Kf/Km$ set from 10^{6} to 10^{5}. When the ratio ($Kf/Km$) is smaller than 10^{4}, the total flow will be increased obviously. For the value of $Kf/Km$ changing from 10^{5} to 10^{6}, the matrix cells in the continuum model may not contribute too much to the total flow through the fractured media. However, when the ratio ($Kf/Km$) is smaller than 10^{5}, the matrix cells would participate actively in flow in simulation and bear considerable part of flow movement.

At the same time, the value of $Km$ is fixed to $1.0\xd710\u22127\u2009m/s$ to analyze the influence of different $Kf$, corresponding to the ratio ($Kf/Km$) varying from 10^{2} to 10^{6}, on the contribution of matrix cells to the total flow. As is shown in Figure 8(b), with the ratio (*K _{f}*/

*K*) changing from 10

_{m}^{4}to 10

^{6}, the slope of total flow growth is almost a constant, implying that the flow in fractures is able to represent the total flow in fracture media including fractures and matrices. When the value of $Kf/Km$ is less than 10

^{4}, the relationship between $Qtotal$ and $Kf/Km$ is deviated from the linearity due to the participation of matrix cells in the total flow. Combining the simulation results under two settings, it is found the conductivity of fracture cells should keep 5 orders of magnitude more than that of matrix cells to ignore the contribution of matrix cells in total flow by using MODFLOW code to simulate flow in 3D FCM models. When the ratio ($Kf/Km$) decreases to less than 10

^{5}, matrix cells would actively participate in the total flowing through the fracture network.

As mentioned above, due to the truncation of isolated fractures, the fracture volume occupied ratio obtained from different realizations of fracture networks in the same case is variable. As a consequence, correlation analysis was applied to quantify the relationship between the fracture volume occupied ratio and the value of $Ke$ for comparison of case 1 with other cases. Note that the both conductivities of $Kf=1.0\xd710\u22122\u2009m/s$ and $Km=1.0\xd710\u22127\u2009m/s$ are set in the flow simulation, representing a typical case of $Kf/Km$ assigned to be 10^{5}.

To examine the stability of realizations generated in the same case, the statistics of $Ke$ in each case were listed in Table 2. It can be found that cases 1 and 2 and cases 4 and 5 have an approximate coefficient of variation (CV). Comparatively, only the CV for case 3 is the smallest which apparently decreased from 22% or so for cases 1-2 and 4-5 to 8.19% for case 3. The results imply that the fracture density plays a significant role for the stability of the method proposed in this study for generating fracture networks. As the fracture volume occupied ratio or the fracture density increases, the backbone of the fracture network is increasingly approaching the original fracture network generated, and then, the value of $Ke$ for each realization tends to be stable.

The correlations between the fracture volume occupied ratio ($\delta $) and the equivalent hydraulic conductivity ($Ke$) of fracture networks obtained from the 200 realizations of generated fracture networks for cases 1-5 are shown in Figure 9, respectively. The coefficients of determination ($R2$) are all greater than 0.825, indicating that there is a good overall positive correlation between $Ke$ versus $\delta $ in all 5 cases. Compared with case 1, cases 2 and 5 have greater $Ke$ due to the existence of the larger length of fractures. Thus, the fracture length is not only a key parameter in fracture connectivity but also an important feature in fracture conductivity. Case 3 is similar to case 1 except for the fracture volume occupied ratio. A comparison of case 3 and case 1 shown in Figure 9(b) indicates that the realizations of fracture networks with higher fracture density reflect more consistency in $Ke$. It can be concluded from the best straight lines for two cases that the relationships between the fracture volume occupied ratio and the equivalent fracture network conductivity should be nonlinear.

Figure 9(c) shows that the fracture orientation affects the hydraulic conductivity. Case 4 has the lowest $Ke$ and $R2$ among all 5 cases. Although the wider distribution range of orientation has advantage in the intersection of fractures, the fracture volume occupied ratio extending direction is not beneficial for the flow movement ability of the whole system. It seems that the $Ke$ is increased when the orientation of fractures becomes closer to the general flow direction. The slope of fitted linear curves for cases 2 and 3 is obviously larger than that for the other 3 cases in Figure 9. This indicates that increasing the average fracture length (e.g., case 2) and the number of fracture grids along flow direction (e.g., case 3) can lead to apparent improvement in equivalent hydraulic conductivity of the whole fracture media. The length and density of fractures are the two important parameters associated with both connectivity and conductivity of the overall fracture networks. It is noteworthy that there is no connectivity of fractures at all if the fracture length decreases from 50 m to 20 m for case 1; this also points to the importance of the fracture length to flow in fracture networks.

To quantify the permeability of fractures set to follow the cubic law in our simulations, the hydraulic conductivity of fractures is only decided by fracture aperture. Therefore, in terms of different scales, the mapping accuracy has the most significant influence on the simulation results. In our simulations, we assume a domain with dimensions $200\u2009m\xd7100\u2009m\xd750\u2009m$, and the resolution of grids is $1\u2009m\xd71\u2009m\xd71\u2009m$. This grid resolution setting may magnify the connectivity in fracture networks when we try to map micro and narrow fractures in a large-scale field (i.e., fractures with a centimeter-level aperture), leading to inaccurately characterizing fracture network permeability. In 2D, it is possible to improve the simulation accuracy by improving the grid resolution, while in 3D, improving the grid accuracy will significantly result in orders of magnitude more computational cost. However, it is worth mentioning that kilometer-scale flow and transport in 3D have been simulated accurately by using FCM method with a larger grid size [30, 40].

## 4. Summary and Conclusions

In this study, a fractured continuum model (FCM) approach was developed for generating fracture networks according to the geometry characteristics, including the distributions of fracture location, density, orientation, length, height, and width. The main results of this paper can be summarized as follows:

- (i)
The generated fracture networks were converted into 3D continuum grids for the simulation of flow in fractured media. Simultaneously, the isolated and dead-end fractures are truncated and eliminated from the originally generated fracture networks. Length, orientation, and density of fractures all play important roles in the connectivity of fractures to intersect each other in the domain

- (ii)
Traditional hydrologic simulator MODFLOW was applied to simulate the flow in the fracture network when the generated fracture network is depicted in continuum grids. For a fracture network, $Kf$ and $Km$ are defined for the fracture and matrix conductivities, respectively. Then, the proper ratio between $Kf$ and $Km$ is achieved at approximately 10

^{5}. It is found that the contribution of fluid flow in matrix can be ignored when $Kf$ is 5 orders of magnitude more than $Km$ in fracture media - (iii)
Five different cases were used to analyze the equivalent hydraulic conductivity ($Ke$) of fracture networks with different generation parameters with the ratio between $Kf$ and $Km$ fixed in 10

^{5}. It is shown that fracture density plays a key role in the stabilization of $Ke$ for the generated fracture networks. Analysis of generated fracture networks shows the highest fracture density can lead to the smallest coefficient of variation (CV) of multirealizations and the highest coefficient of determination ($R2$) of fracture density vs. $Ke$. Furthermore, the value of $Ke$ increases with the enlarging length of fractures and closing orientation of fractures to the general flow direction

It has to be admitted that the FCM approach is a still-evolving solution for solving fracture network flow and transport problems. Like other emerging approaches, it has pros and cons. Its advantage is that the approach combines the main features of traditional DFN and SC methods. The approach can solve the hydraulic exchange between the fracture and the matrix while preserving the geometric characteristics of the fracture network. Its disadvantage comes from the connectivity between the fracture and the matrix in the grids. The enhanced connectivity changes the total flow of the fracture network, which will be overestimated. Therefore, the significance of the research on the ratio of hydraulic conductivity is that the contribution of the matrix to the total flow can be accurately controlled in the further work. For the FCM-based fluid flow and solute modeling in real-world applications such as those in Figure 1, multidisciplinary and multiscale approaches including acoustic, seismic, electric resistivity tomography, and core data [41, 42] are needed to jointly investigate the heterogeneity and its overall permeability of a fractured aquifer composed of fracture and matrix in future studies.

Another limitation of this study is that the effect of grid size on the simulation results is not discussed. The finer grid size setting will lead to more accurate simulation results. However, it will be found that the approach becomes computationally intensive because of massively increased grids. This is a topic worthy of further discussion. The general conclusion of this study is that the approach shows the possibility of applying the commonly used finite difference method in porous media to simulate the flow in 3D fracture networks. It provides a stepping stone for further research on fluid flow and solute transport problems in more complex fracture networks, especially for the field-scale real-world applications.

## Nomenclature

- $A$:
Subregion of $R$ in homogenous Poisson process

- $An,Ak,Ai$:
Disjoint subregions of

*R*in homogenous Poisson process- CV:
Coefficient of variation

- $J$:
Hydraulic gradient

- $Ke$:
Equivalent hydraulic conductivity (m/s)

- $Kf$:
Fracture conductivity (m/s)

- $Km$:
Matrix conductivity (m/s)

- $Kxx$:
Hydraulic conductivity tensor in $x$ direction (m/s)

- $Kyy$:
Hydraulic conductivity tensor in $y$ direction (m/s)

- $Kzz$:
Hydraulic conductivity tensor in $z$ direction (m/s)

- $lmean$:
Mean of fracture length (m)

- $N$:
Number of points

- $qs$:
Sink/source term (m

^{3}/s)- $Q$:
Total grid flow moving from upstream to downstream boundaries (m

^{3}/s)- $Qtotal$:
Total flow discharge moving from upstream to downstream boundaries in simulation (m

^{3}/s)- $R$:
Region of interest in homogenous Poisson process

- $R2$:
Coefficient of determination

- $Uj$:
Random variable in homogenous Poisson process.

## Data Availability

The data can be available by contacting the corresponding author or the first author.

## Conflicts of Interest

The authors declare that they have no conflict of interest.

## Acknowledgments

We gratefully acknowledge the financial support from the National Key Research and Development Plan of China (2019YFC1805302) and the National Natural Science Foundation of China (U2167212 and 42272277). We are also grateful to the High Performance Computing Center of Nanjing University for providing the computing facility. Especially, the authors are profoundly grateful to the editor Prof. Jean Borgomano and two anonymous reviewers, whose constructive comments and conscientious suggestions led to significant improvement of the manuscript.