Effective Permeability in Fractured Reservoirs: Discrete Fracture Matrix Simulations and Percolation‐Based Effective‐Medium Approximation

Fractured reservoirs are complex and multi‐scale systems composed of matrix and fractures. Accordingly, modeling flow in such geological media has been a great challenge. In this study, we investigated the effect of scale as well as matrix and fracture network characteristics on the effective permeability (keff) in matrix‐fracture systems under fully saturated conditions. We generated fracture networks, embedded within a matrix of permeability of 10−18 m2, with fracture lengths followed a truncated power‐law distribution with exponent α = 1.5, 2.0, and 2.5. We set fracture permeability equal to 10−16, 10−14, and 10−12 m2 and numerically simulated fluid flow to determine the keff at six fracture densities for 36 fractured reservoirs. Results showed that the effect of α and scale on the keff became more significant as the contrast between matrix and fracture permeabilities increased. We also fit the percolation‐based effective‐medium approximation (P‐EMA) to the simulations and optimized its two parameters critical fracture density and scaling exponent. Results exhibited that both P‐EMA model parameters were scale‐dependent. Through linear regression analysis, we found that the critical fracture density and scaling exponent were highly correlated to other matrix‐fracture system properties and proposed two regression‐based models evaluated using a new six sets of simulations. Comparing the estimated keff values with the simulated ones demonstrated the reliability and predictability of the P‐EMA. The matrix‐fracture systems studied here were finite in size. We also showed that one may extend results to infinitely large reservoirs using the P‐EMA framework.

fracture network (DFN) approach was widely employed to model flow and transport (Berrone et al., 2018;Dverstorp et al., 1992;Medici et al., 2021;Medici & West, 2022;Painter & Cvetkovic, 2005).In the DFN method, governing equations of flow or transport are numerically solved in individual fractures.Although ingenious, this approach overlooks the permeability contribution from the rock matrix, assuming it to be negligible.Although incorporating matrix flow is computationally demanding, some researchers (e.g., Deutsch, 1989;Pierce et al., 2018) studied flow and transport in matrix-fracture systems.For example, Sweeney et al. (2020) introduced the upscaled discrete fracture matrix (UDFM) model, an alternative to the traditional DFN approach, to capture interactions between fractures and surrounding rock matrix more accurately.Those authors demonstrated their model applicability to complex heterogeneous fractured media and validated it by comparing with numerical and analytical benchmarks.
In addition to numerical methods, theoretical models were developed to study k eff in fracture networks and matrix-fracture systems.One of the early models is the power-law averaging (Deutsch, 1989;Journel et al., 1986).For instance, Zanon et al. (2002) applied power-law averaging to formations constructed of high-permeability sandstone and low-permeability shale grids.Those authors found that the power-law exponent (ω) depended on geological properties, such as the percentage of sandstone and anisotropy ratio.Years later, de Dreuzy et al. (2001b) applied the power-law averaging to permeability in two-dimensional fracture networks with apertures lognormally distributed.They reported that the exponent ω varied with network size, fracture length powerlaw distribution exponent, and fracture density.In another study, Mourzenko et al. (2011) proposed a heuristic model based on the asymptotic behavior of permeability in three-dimensional fracture networks, consistent with numerical data simulated over a wide range of network densities.However, their model is not applicable to matrix-fracture systems where matrix contributes to fluid flow at very low fracture densities (ρ < ρ c where ρ and ρ c are respectively fracture density and its critical value).More recently, following the work of Saevik et al. (2013), Ebigbo et al. (2016) evaluated various effective medium-based models, such as symmetric and asymmetric self-consistent, differential, and Maxwell, as well as the heuristic model of Mourzenko et al. (2011).They numerically simulated effective permeability in three-dimensional fractured rock masses composed of spheroidal fractures with different aspect ratios.In their simulations, Ebigbo et al. (2016) considered both mono-and poly-disperse fracture networks and varied the matrix permeability to fracture permeability ratio from 1.2 × 10 −7 to 4.8 × 10 −5 .They found good agreement between theoretical estimations and numerical simulations for the self-consistent effective-medium approximation.Ebigbo et al. (2016) also reported that the heuristic model of Mourzenko et al. (2011) was accurate, particularly for mono-disperse networks.
Percolation-based effective-medium approximation (P-EMA) proposed in physics literature by McLachlan (1987McLachlan ( , 1988) ) provides another theoretical framework to study fluid flow in two-component systems.However, to the best of our knowledge, neither has it been applied to matrix-fracture systems nor has it been attempted to model effective permeability in fractured reservoirs.Therefore, the main objectives of this study are to: (a) apply the P-EMA to fit k eff as a function of ρ, (b) investigate how the P-EMA parameters (critical fracture density ρ c and scaling exponent t) vary with fractures and reservoir properties, (c) study effects of heterogeneity on the k eff in the matrix-fracture systems, and (d) address the effect of finite size on fractured reservoir and simulation of fluid flow through them.

Numerical Simulations
For the DFN simulations, we applied a parallelized computational platform called dfnWorks developed at the Los Alamos National Laboratory (Hyman et al., 2015).We utilized the Upscaled Discrete Fracture-Matrix (UDFM) model developed by Sweeney et al. (2020), which incorporates the effect of rock matrix and its contribution to the k eff .In the following sections, we explain the fracture networks and flow simulations in further detail.

Fracture Network Generation
Within the dfnWorks computational platform, fracture networks are generated using the dfnGen package, which utilizes two different libraries: (a) feature rejection algorithm for meshing (FRAM) and (b) Los Alamos grid toolbox (LaGriT).Each DFN is constructed so that all features in the network, for example, length of intersections between fractures and distance between lines of intersection of a fracture are greater than a user-defined minimum length scale (Hyman et al., 2015).Like most DFN modeling techniques, one of the key advantages of the Observational data from fractured media in nature show that the probability density function of fracture lengths (l) is broad and typically conform to the following truncated power-law probability density function (Bonnet et al., 2001;de Dreuzy et al., 2012;Hyman et al., 2018): where l min and l max are the minimum and maximum fracture lengths, respectively, c l is a constant coefficient (a normalization factor), and α is the exponent that controls the frequency of fracture length and, therefore, the heterogeneity of the network (de Dreuzy et al., 2012).The value of α varies typically between 1 and 3 (Bonnet et al., 2001).
Within the dfnWorks framework, each fracture has a specific orientation, which is sampled from the Fisher distribution where β is the dip angle (the mean orientation vector), κ is a concentration parameter controlling the uniformity and heterogeneity of the fracture orientation.The range of κ is wide and typically varies from 0 to ∞. κ values closer to zero tend to generate a more uniform distribution with all fracture orientations equally likely to be generated, while larger κ values tend to generate fracture orientation clustered around the mean orientation.
Following Vermilye and Scholz (1995), we set l min = 0.02 m and l max = 20 m, in accord with the experimental observations reported for Florence Lake.Three different α values were considered (α = 1.5, 2.0, and 2.5) to cover the experimental range reported by Bonnet et al. (2001).We also assumed that fracture length and aperture data collected from Florence Lake were correlated (Vermilye & Scholz, 1995), also highlighted by Olson (2003).To further generate a polydisperse fracture network that matched the experimental data, we set κ = 0.1 and β = 45 °.
Although seismic data were used to extract and characterize fracture networks (Thachaparambil, 2015), capturing all properties of matrix-fracture systems in subsurface is extremely difficult.Therefore, our aim was not generating fracture networks observed exactly in nature but statistically representative ones.

Meshing
We applied an octree-refined continuum mesh (Sweeney et al., 2020) to adequately capture the geometrical features of the generated matrix-fracture systems.The generated matrix-fracture systems were mapped onto a uniformly discretized hexahedral mesh to account for the presence of rock matrix which serve as the background of the mesh.After mapping the DFN onto the hexahedral mesh, the mesh was then binarized into fracture and matrix cells with cells intersected by fracture tagged "fracture cells" and cells not intersected by fracture tagged "matrix cells."The resolution of the mesh depends on the proximity to individual fractures within the fracture network.Higher mesh resolutions were employed near the fractures in the system which enabled us to capture the dynamic processes that occurs between the fractures and the surrounding matrix.The octree method utilizes two user-defined parameters: (a) edge length of the original hexahedral mesh before refinement (l) and (b) the number of refinement level (r) in the final octree mesh.In our simulation setup, l was set equal to L/10 and r equal to 2. Following Sweeney et al. (2020), who reported that the value of r is less likely to significantly impact k eff , we set r = 2 to reduce the computational cost involved in meshing while ensuring accuracy in k eff calculation.

Fluid Flow
Following the binarization of the system into fracture and matrix cells, their respective permeabilities (i.e., k f and k m ) were upscaled into their respective cells.The upscaled properties were incorporated in the fluid flow simulation, and accordingly the value of k eff was determined at different ρ values between 0 and 1. Fluid flow was simulated using PFLOTRAN, a parallel subsurface flow and reactive transport finite volume code (Lichtner et al., 2015), under fully saturated conditions.The PFLOTRAN employs Richards' equation, extensively applied in subsurface hydrology for modeling flow under single phase, variably saturated, isothermal, and steady stay conditions, to compute the outlet volumetric flow rate.Pressure boundary conditions at the inlet and outlet was set at 1.1 × 10 6 and 1.0 × 10 6 Pa, respectively.We then numerically invert Darcy's equation to calculate the k eff of individual DFN in the network: where q is the Darcy flux, calculated by dividing the volumetric flow rate Q, which is computed by PFLOTRAN, by the area of the system (L × L), ∆p is the change in pressure from the inflow to outflow boundary, and L is the system size.
In this study, the k m was set equal to 10 −18 m 2 in accord with shale matrix permeability values reported by Best and Katsube (1995) and Wang et al. (2009).Three fracture permeability values k f = 10 −16 , 10 −14 , and 10 −12 m 2 corresponding to log(k f /k m ) = 2, 4, and 6 were considered.These values, k m and k f , were upscaled in their appropriate cells, that is, matrix or fracture cell.
For a given complete simulation, we extracted the value of k eff and computed the corresponding fracture density (ρ), that is, the fraction of volume occupied by fracture within the matrix-fracture system to the overall volume of the system, which were subjected to further analyses.The equation for fracture density is given by: where v f is the volume occupied by fracture in the system while v t is the total volume of the system.
As stated earlier, we conducted iterations of the matrix-fracture system generation and fluid flow simulations to obtain statistically reliable results.We started with 10 iterations at each ρ, and then computed the average, variance and 95% confidence intervals (CIs) of the k eff .The average k eff was then plotted against the corresponding mean ρ with the 95% CIs.Each individual k eff simulation was next added to the plot, and those fallen outside the 95% CIs were eliminated from further analyses.After that, we calculated the standard deviation of the iterations that fell within the 95% CIs.More iterations were conducted, and the above process was repeated until the standard deviation no longer changed with the number of realizations.

Percolation-Based Effective-Medium Approximation
Percolation-based effective-medium approximation (P-EMA) is an upscaling technique from statistical physics originally proposed by McLachlan (1987McLachlan ( , 1988) ) for media with low-and high-conductivity components.It includes percolation theory and effective-medium approximation as its special cases (Ghanbarian & Daigle, 2016).Although previously applied to study electrical conductivity, thermal conductivity, permeability and Young's modulus in binary composites (Deprez et al., 1988;McLachlan, 2021), its applications to porous rocks and reservoirs have been very limited.To the best of our knowledge, the P-EMA has not yet been applied to model the k eff in matrix-fracture systems.
Within the P-EMA framework, the relationship between   and   is given by where ρ c is the critical fracture density and t is the scaling exponent.In Equation 5, k eff is implicitly explained in terms of ρ.Note that when  5. The P-EMA approach is rather general and includes several other models as its special cases (Ghanbarian & Daigle, 2016;Sadeghi et al., 2018).For instance, Equation 5 reduces to the weighted arithmetic mean when ρ c = 0, to the weighted harmonic mean when ρ c = 1, to the geometric mean when ρ c = ρ = 0.5, and to the Bruggeman (1935) model when ρ c = 1/3.Equation 5 also reduces to the power-law scaling relationship   ∝ ( − )  from percolation theory when the matrix permeability is zero (k m = 0).
Rearranging Equation 5gives ρ explicitly as a function of k eff as follows Note that for ρ < ρ c , k eff is mainly controlled by the rock matrix, while for ρ > ρ c , k eff is dominated by the fracture network and its properties.In this study, we fit Equation 6to the averaged ρ−k eff curves using the nonlinear least square optimization method, optimized the P-EMA model parameters that is, ρ c and t, and explored the relationship between the matrix-fracture system properties and the P-EMA model parameters.
The P-EMA approach was experimentally validated at the lab scale in the literature (Deprez et al., 1988;McLachlan, 1986McLachlan, , 1988McLachlan, , 1996)).For instance, McLachlan (1988) fit the P-EMA model to two-dimensional electrical conductivity data experimentally measured on a plate consisted of overlapping insulating squares (continuum percolation) and found a percolation threshold of 0.39.This value is very close to 0.407 reported by Dubson and Garland (1985) who measured the electrical conductivity of a sheet of aluminized Mylar in which insulating squares were inscribed.Deprez et al. (1988) also demonstrated applications of the P-EMA model to electrical and thermal conductivities as well as Young's module for a series of sintered nickel samples.More recently, Ghanbarian and Daigle (2016) used the P-EMA model to characterize the thermal conductivity of porous media under fully and partially saturated conditions and found very good fit to experimental data.

Multiple-Linear Regression Analysis
Multiple-linear regression analysis has been widely used to present simple and statistical relationships between a target and some input variables.To make the P-EMA model predictive, we applied stepwise multiple-linear regression analysis using the statistical software SAS and linked the P-EMA parameters t and ρ c to other matrix-fracture system characteristics.To establish regression-based models, some new features from existing ones were engineered, and independent variables α, L,

and 𝐴𝐴
√ ∕ were used.The Jarque-Bera test was applied to examine whether or not the distribution of input variables was normal.For training the models, all the data including 36 simulations (summarized in Table 1) were utilized.After training, the normality of the error term was evaluated using the Quantile-Quantile plot.
To assess the reliability and predictivity of the developed regression-based models, we conducted new six simulation sets based on outcrop data of three naturally fractured sites that is, Culpeper quarry (Vermilye & Scholz, 1995), Hornelen1 bed (Azizmohammadi & Matthäi, 2017), and Kilve bed (Lei et al., 2015).In the first four simulations, we set l min = 0.5 m, l max = 13 m, following Vermilye and Scholz (1995), k m = 10 −18 m 2 , and L = 100 m.We also used α = 2.25, 2.25, 1.75 and 1.75 and log(k f /k m ) = 5, 7, 5, and 7 in simulation sets 1, 2, 3, and 4, respectively.In the simulation set 5, we assumed l min = 1.1 m, l max = 10 m, L = 18 m and α = 2.83, following Azizmohammadi and Matthäi (2017), k m = 10 −15 m 2 , and log(k f /k m ) = 9.In the simulation set 6, we set l min = 0.3 m, l max = 3 m, L = 6 m and α = 2.37, following Note.k m = 10 −18 m 2 consistent with shale matrix permeability.L = 22.5, 30, 50, and 70 m were selected so the system size to be greater than the maximum fracture length l max = 20 m and to avoid one single long fracture cross the entire system.

Table 1
The Optimized Values of the P-EMA Model Parameters for Different Matrix-Fracture Systems Generated in This Study Lei et al. (2015), k m = 10 −15 m 2 , and log(k f /k m ) = 9.For all simulation sets, we assumed κ = 0.1 and β = 45° in Equation 2.

Results and Discussion
In this section, we present the results of the simulated k eff as a function of ρ for different matrix-fracture systems, fitting the P-EMA model to the k eff -ρ simulations for 36 matrix-fracture systems, relationship between the P-EMA model parameters and matrix-fracture system properties, and the extension of our results to an infinitely large fractured reservoirs.
We demonstrate several DFNs generated in this study in Figure S1 in Supporting Information S1.As can be seen in Figure S1 in Supporting Information S1, the networks generated with α = 1.5 tended to have longer fractures and better connected compared to the networks with α = 2 and 2.5.As Figure S1 in Supporting Information S1 shows, the DFN generated with α = 2.5 composed of shorter fractures, which ultimately resulted in relatively high percolation threshold because more fractures are needed to form a sample spanning cluster percolating the network.
To minimize uncertainties in the stochastically generated matrix-fracture systems and to obtain representative and statistically reliable results, we had multiple realizations.The number of realizations varied from one system to another depending on the value of α and the system size (L).The total number of iterations required to achieve statistically reliable results for each matrix-fracture system realization depended on the level of heterogeneity, which is largely controlled by α and L. For α values of 2.5 and 2.0, a minimum of 10 iterations were conducted for each ρ in a complete set of simulations comprising 6 data points (six simulated permeabilities against six fracture densities on the k eff −ρ curve).However, due to the increased heterogeneity near the critical fracture density, which is the minimum ρ required for percolation to occur within the system, at least 20 iterations were conducted at intermediate fracture densities where the transition of k m to k f happens.Consequently, a single matrix-fracture system realization generated at α values of 2.5 and 2.0 and L ≥ 50 m had a total of at least 80 iterations.For smaller values of L, where we observed an increased heterogeneity, such as L = 22.5 and 30 m, we conducted at least 20 iterations for each data point irrespective of the value of α.The DFNs generated with α = 1.5, as discussed earlier, have a higher tendency to generate longer fractures in the network, which significantly increased the level of heterogeneity in the matrix-fracture system relative to the DFNs generated with α = 2.5 and 2.0.Therefore, for α = 1.5, we had at least 20 iterations for each ρ value in a complete set of simulation, and a minimum of 30 iterations were conducted around the critical fracture density.In total, we conducted over 5,000 simulations for different matrix-fracture systems summarized in Table S1 in Supporting Information.

Effect of α on k eff in Matrix-Fracture Systems
Figure 1 shows the average k eff against the ρ for different α values.At lower ρ values where the hydraulic properties of the matrix-fracture system are controlled by the rock matrix, k eff remains nearly constant.At some intermediate fracture density, however, fluid finds a conductive pathway through the fracture network, which results in a significant increase in k eff .At higher ρ values, the k eff −ρ curve becomes flattened toward k f yielding a sigmoidal shape for log(k f /k m ) > 2 (Figure 1).Note that fracture densities as small as 0.1 (Karaman & Carpenter, 1997) and as large as 0.6 (Clarke & Burbank, 2011) have been reported in the literature.
Results presented in Figure S1 in Supporting Information S1 confirm that lower α values rendered longer fractures within the fracture network, which significantly impacted connectivity and consequently the k eff −ρ curve.This is consistent with the results of Berkowitz et al. (2000) and others who reported that the connectivity of fracture networks was dependent on the exponent α.
We identified two distinct regions: (a) the matrix-controlled region where ρ < ρ c and (b) fracture-controlled region where ρ > ρ c . Figure 1 shows that for α = 1.5 the transition from the matrix-controlled region to the fracture-dominated one occurred at a lower fracture density ρ compared to α = 2.5.Additionally, as observed in Figure 1, the transition from k m to k f becomes more significant as the log(k f /k m ) value becomes greater.Zhu et al. (2021) investigated how geometrical properties of fracture networks, such as length, orientation, aperture, and position of fracture centers, affect macro-scale flow properties in shale-like formations.They found that as the fracture length power-law distribution exponent (α) increased, the global efficiency of the network decreased, that is, the connectivity of the system decreased with increasing values of α, and more fractures were required to achieve percolation.They also reported that flow rate increased as α decreased, which is consistent with our results shown in Figure 1.In their study, however, Zhu et al. (2021) neither report effective permeability nor address the effect of scale on flow rate.They also did not provide a theoretical framework to explain how flow rate varies with other matrix-fracture system characteristics.

Effect of log(k f /k m ) on k eff in Matrix-Fracture Systems
In Figure 2, we show the value of k eff , averaged over several realizations, as a function of ρ for log(k f /k m ) = 2, 4, and 6 and α = 1.5, 2, and 2.5.Our results clearly indicate that the transition from the matrix-controlled region (ρ → 0) to the fracture-dominated one (ρ → 1) was more significant with increasing the value of log(k f /k m ), as illustrated in Figure 2.This means that in a connected fractured medium, the magnitude at which fractures control flow depends on log(k f /k m ).Comparing Figures 2a-2c also shows that such a transition happened at a larger fracture density as the value of α increased from 1.5 (Figure 2a) to 2.5 (Figure 2c).Our results agree with the results of Hyman et al. (2018) who reported that the magnitude of effective permeability increase around the percolation threshold was greater in matrix-fracture systems of greater log(k f /k m ).Results presented in Figure 2 are also consistent with those reported by Ebigbo et al. (2016) who found that in systems with small perturbations in matrix and fracture permeability values, the simulated k eff −ρ data exhibited nearly a linear trend with smooth transition from k m to k f .However, the percolation threshold was more distinctive as perturbations increased.

Effect of Scale L on k eff in Matrix-Fracture Systems
Understanding the scale dependence of intrinsic properties of porous media, such as permeability, is necessary, particularly for upscaling purposes.For a given simulation with consistent input parameters at different values of L, at ρ < ρ c we observed that k eff remained almost the same for L = 22.5, 30, 50, and 70 m.At ρ > ρ c , we found that lower L values exhibited higher k eff and k eff tended to decrease as L increased, which is consistent with the findings of Lei et al. (2015) who investigated the effect of scale in two dimensions.de Dreuzy et al. ( 2001a) also reported similar behavior for two-dimensional DFNs generated with α < 3, consistent with the range of α in this study.From our analyses, we observed that percolation occurred relatively early at low L values, and the percolation threshold tended to increase with increasing L, in accord with the finite-size scaling analysis (Stauffer & Aharony, 2018) as we discuss in further detail later.We further discuss the scale dependence of k eff for both finite and infinitely large systems in Section 5.6 where we extrapolate our results to infinitely large reservoirs.

Percolation-Based Effective-Medium Approximation
Using the Curve Fitting Toolbox of MATLAB, the P-EMA, Equation 6, was fit to the ρ−k eff simulations averaged over various realizations.The optimized values of the P-EMA parameters, ρ c and t, for the matrix-fracture systems studied here are summarized in Table 1.As an example, the P-EMA fits to the data corresponding to L = 50 m and different values of log(k f /k m ) and α are shown in Figure 3.The high values of R 2 (=0.99) reported in Table 1 imply that the P-EMA (Equation 6) fit the data well.Using the data reported in Table 1, we plotted the critical fracture density ρ c against the exponent α for each system size and observed that the smaller the α value the lower the ρ c , meaning that the lower α values resulted in early percolation in the fracture network.We found the lowest value of ρ c for the matrix-fracture systems with α = 1.5, followed by α = 2.0, while the matrix-fracture systems with α = 2.5 had the highest value of ρ c .Our results, however, are not consistent with those of de Dreuzy et al. (2000) who reported that by increasing α the percolation threshold in three-dimensional fracture networks decreased.de Dreuzy et al. (2000) argued that such a trend may be attributed to truncation effects because in their study elliptical fractures, truncated by system sides, had an internal characteristic length less than the original one.In addition to that, in their study α ranged between 2.5 and 5, while in our study between 1.5 and 2.5 in accord with the experimental range reported for naturally fractured sites (Bonnet et al., 2001).Recall that α controls the frequency of fracture length; the lower the α value, the greater the number of longer fractures (Figure S1 in Supporting Information S1).Therefore, the correlation between ρ c and α seems reasonable.
The scaling exponent, t, in Equations 5 and 6 was observed to decrease as the value of α increased.For example, the average value of t for all matrix-fracture systems with α = 1.5 is 2.07 but decreases to 2.02 and 1.85 for α = 2.0 and 2.5, respectively.Additionally, increasing the system size, L, resulted in a decrease in the t value.For L = 22.5, 30, 50, and 70 m, the corresponding average values of t are 2.15, 2.05, 1.87, and 1.84, respectively.Moreover, our analyses revealed that t exhibits an inverse relationship with log(k f /k m ).The average values of t for log(k f /k m ) = 2, 4, and 6 are 2.19, 1.92, and 1.83, respectively.Since log(k f /k m ) controls the shape of the k eff −ρ curve and it is inversely correlated to t, one may expect the lower value of the scaling exponent t to result in sharper increase in the k eff .

Regression-Based Models for Estimating k eff
In this section, we further investigated the relationship between the P-EMA model parameters (ρ c and t) and other matrix-fracture system properties through stepwise multiple-linear regression analysis.We found that Regression-based results showed that ρ c was significantly and statistically linked to α and L (p-value of <0.0001) consistent with our previous results stated earlier.Although positive correlation between ρ c and α was also reported by Mourzenko et al. (2005) for fracture networks with α < 3, Sahimi and Mukhopadhyay (1996) found an inverse relationship.They, however, studied the scale dependence of percolation threshold in networks with long-range correlations.Drawing upon our previous analysis of the impact of L on k eff , we found that when holding all other input parameters constant, an increase in L led to an increase in the percolation threshold.This implies that more fractures are required in larger systems to form a sample spanning cluster compared to smaller ones.This supports the results obtained from our regression analysis, which confirms that L has a significant influence of ρ c .
Our regression-based results showed that the independent variables α, log(k f /k m ), and L statistically and significantly contributed to the dependent variable t (p-value < 0.0001).We also found that t was negatively correlated to α, L, and log(k f /k m ) (see Equation 8).Such dependencies suggest that t is influenced by geometrical properties of the matrix-fracture system as well as the hydraulic properties of the system.
The accuracy of Equations 7 and 8 and the estimated k eff at different ρ values are presented in Figure 4 for the six simulation sets.We first estimated values of ρ c and t via Equations 7 and 8, respectively.Equation 6 was then used  2.
to estimate the k eff at various ρ values.The root mean square log-transformed error (RMSLE), relative absolute error (RAE), and correlation coefficient (R 2 ) values are also reported in Figure 4.As can be seen, the estimations matched the simulations well with 0.11 < RMSLE < 1.96, 16.5% < RAE < 75.3%, and 0.98 < R 2 < 0.99, which confirmed the predictability of the P-EMA model.

Extrapolation to Infinitely Large Fractured Reservoirs
Modeling fluid flow within the actual size of a matrix-fracture system in nature, which may be practically considered as an infinitely large medium, would be challenging due to computational costs.In this section, we apply concepts of finite-size scaling analysis (Stauffer & Aharony, 2018), to approximate the value of ρ c in fractured reservoirs whose dimensions are infinitely large.We also use an empirical exponential relationship proposed by Matyka et al. (2008) to estimate the value of t in infinitely large matrix-fracture systems.The predicted values of ρ c and t reported were obtained from Equations 7 and 8 respectively.For all simulations we assumed κ = 0.1 and β = 45°.

ρ c for Infinitely Large Systems
Within the finite-size scaling analysis framework, one can explain the scale dependence of ρ c as follows (Stauffer & Aharony, 2018): where ρ c (L) and ρ c (L → ∞) are the critical fracture densities for a finite-and infinite-sized matrix-fracture systems, respectively, C is a constant coefficient whose units depends on the system units, and ν is the correlation length exponent equal to 0.88 in three dimensions (Hunt et al., 2014).Finite-size scaling analysis has been widely applied in the literature (Sahimi, 2011).For example, Ji et al. (2004) let C and ν to be fitting parameters, fit Equation 9 to simulations on two-dimensional fracture networks and reported C = 0.15 and ν = 1.27.The latter is slightly less than the universal ν value in two dimensions that is, 1.33 (Hunt et al., 2014).Mourzenko et al. (2011) also applied a model similar to Equation 9 to extrapolate permeability to infinitely large networks. and the fit of Equation 9 to the corresponding data.For fitting purposes, we used the Curve Fitting Toolbox of MATLAB.As can be seen, Equation 9 fit the data well with an average R 2 = 0.98.We found that the optimized value of C, reported in Figure 5a and Figure S2 in Supporting Information S1, decreased with increasing log(k f /k m ).For instance, for log(k f /k m ) = 2, 4, and 6, C was respectively 6.93, 6.58, and 6.5 when α = 1.5, while 4.60, 3.68, and 3.20 when α = 2.5.
The optimized value of ρ c (L → ∞) for all nine matrix-fracture systems is presented in Table 2. Similar to our previous results, we found that the value of ρ c (L → ∞) increased with the increase of α. Results tabulated in Table 2 also show that the value of ρ c (L → ∞) in the systems with the same α (particularly α = 1.5 and 2) did not change from log(k f /k m ) = 2 to 6, which means that ρ c (L → ∞) depends more on the value of α than log(k f /k m ).
This finding is consistent with our regression analysis results; no statistical evidence of a relationship between ρ c and log(k f /k m ).Matyka et al. (2008) proposed an empirical relationship to study the scale dependence of tortuosity and to determine its value for an infinitely large medium.Similarly, we applied the following exponential relationship to explain the scale dependency of the exponent t and extrapolate its value for

t for Infinitely Large Systems
where b and c are two constant coefficients and t inf is the value of t for an infinite-sized matrix-fracture system.In Equation 10, as L approaches infinity, t tends to t inf .
We plotted t against L and applied Equation 10 to fit the data.Results presented in Figure 5b; Figure S3 in Supporting Information S1 show that  7).(b) The scaling exponent t against the system size L for α = 1.5 and log(k f /k m ) = 6.The red line represents the fit of Equation 8to the data.Equation 10 fit the t-L data well with an average R 2 value of 0.99.We also listed the value of t inf for all nine matrix-fracture systems in Table 2.As can be seen, the value of t inf decreased as the value of log(k f /k m ) increased.
We also observed that greater α values corresponded to smaller t inf values, which is in well accord with our regression-based results presented in Equation 8.
The decreasing trend of t with increasing L shown in Figure 5b and Figure S3 in Supporting Information S1 is also consistent with the results of Tremblay and Machta (1989) who theoretically demonstrated that the scaling exponent t is scale dependent and numerically showed that its value decreased as L increased in two and three dimensions.

Estimating the k eff −ρ Curve for Infinitely Large Systems
To extend our results to infinitely large systems, we replaced ρ c and t in Equation 6with the calculated values of ρ c (L → ∞) and t inf reported in Table 2 and determined the k eff at various ρ values via the P-EMA.Results are given in Figure S4 in Supporting Information S1 for α = 1.5, 2.0, and 2.5, and log(k f /k m ) = 2, 4, and 6. Figure S4 in Supporting Information S1 also shows the k eff −ρ curves corresponding to L = 22.5, 30, 50, and 70 m.As can be observed, k eff −ρ curves for L = 22.5, 30, 50, 70 m and infinity and log(k f /k m ) = 2 are almost inseparable (Figures S4a, S4d, and S4g in Supporting Information S1) meaning that the effect of scale on the k eff was negligible for small permeability perturbations.However, due to higher level of heterogeneity, the influence of scale became more substantial as the value of log(k f /k m ) increased from 2 to 6.This means the greater the log(k f /k m ) value, the more profound the effect of scale on k eff in the matrix-fracture systems.We found as L increased, the percolation threshold increased as well.However, for α = 1.5 and log(k f /k m ) > 2 (Figures S4b and S4c in Supporting Information S1), the percolation threshold values for L = 30, 50, and 70 m were similar.
As stated earlier, the networks generated with α = 1.5 had higher frequency of longer fractures (see Figure S1 in Supporting Information S1), which resulted in an earlier percolation compared to those generated with α = 2 and 2.5.
We also found that for α = 1.5 and  log ( ).This is due to the fact that the system size L = 22.5 m is slightly greater than the l max = 20 m in our simulations.For such small systems, there is a higher probability that a single long fracture propagates the entire system, resulting in a smaller percolation threshold.We also observed that the value of log(k f /k m ) impacted the shape of the k eff −ρ curve.For log(k f /k m ) = 2, the k eff kept increasing at ρ ≥ ρ c until it reached k f near ρ = 1 yielded an increasing concave upward trend.However, as the value of log(k f /k m ) increased to 4 and 6, the trend turned into sigmoidal.

Conclusions
In this study, we investigated the effective permeability, k eff , in the matrix-fracture systems under fully saturated conditions by means of numerical simulations and theoretical modeling.To address the effect of scale and matrix-fracture system properties, fluid flow was numerically simulated by solving Richards' equation in 36 fractured reservoirs.To eliminate uncertainties, the simulations were iterated at least 10 times at each fracture density, ρ, with more than 5,000 iterations in total.The simulated k eff −ρ curves were then averaged and fit by the percolation-based effective-medium approximation (P-EMA) and its parameters, ρ c (critical fracture density) and t (scaling exponent), were optimized.We demonstrated that both ρ c and t were scale-dependent.We also found that the effect of scale was more significant in the systems with greater log(k f /k m ) values (=4 and 6).Our numerical simulations indicated that the P-EMA parameters ρ c and t depended on the matrix-fracture characteristics, such as α, log(k f /k m ), maximum fracture length (l max ), and system size (L).Using the stepwise multiple-linear regression analysis, we developed two regression-based models, linked ρ c and t to other matrix-fracture properties and found accurate estimations of k eff via the P-EMA model for six unseen matrix-fracture systems.Our results clearly demonstrated the predictability of the P-EMA.However, the developed regression-based models should be applied with caution to estimate the k eff in matrix-fracture systems substantially different from those studied here.We also extended our simulations to infinitely large fractured reservoirs by determining the values of ρ c and t at L → ∞.Our theoretic approach based on the P-EMA not only provides a framework to estimate k eff in matrix-fracture systems but also incorporates the effect of system size to scale up effective permeability in fractured reservoirs.
is its flexibility in accommodating any statistical survey of a fractured site for fracture network generation.This allows the generation of DFNs that are representative of naturally fractured sites.
which is well met by Equation 5 since under such a circumstance the terms would be zero, equal to the right side of Equation

Figure 1 .
Figure 1.Effective permeability, k eff , versus fracture density, ρ, for (a) log(k f /k m ) = 2, (b) log(k f /k m ) = 4, and (c) log(k f /k m ) = 6 and different α values.In all cases, the system size, L = 50 m and the matrix permeability k m = 10 −18 m 2 .Each data point represents the average over multiple iterations and the error bars correspond to one standard deviation.The impact of α appears to be more significant with increasing log(k f /k m ).

Figure 2 .
Figure 2. Effective permeability, k eff , versus fracture density, ρ, for (a) α = 1.5,(b) α = 2.0, and (c) α = 2.5 and different log(k f /k m ) values.In all cases, the system size, L = 50 m and the matrix permeability k m = 10 −18 m 2 .Each data point represents the average over multiple iterations and the error bars correspond to one standard deviation.

Figure 3 .
Figure 3.The fit of the P-EMA to the k eff −ρ simulations for (a) α = 1.5,(b) α = 2.0, and (c) α = 2.5 and different log(k f /k m ) values.In all cases the system size L = 50 m.The optimized parameters of the P-EMA are summarized in Table2.

Figure
Figure 5a and Figure S2 in Supporting Information S1 show ρ c as a function of   − 1

Figure 5 .
Figure 5. (a) Relationship between critical fracture density (ρ c ) and system size for α = 1.5 and log(k f /k m ) = 6.Fitted red line indicates the finite-size scaling equation (Equation7).(b) The scaling exponent t against the system size L for α = 1.5 and log(k f /k m ) = 6.The red line represents the fit of Equation8to the data.To see the fits to the rest of the cases, see FiguresS2 and S3in Supporting Information S1.
Figures S4b and S4c in Supporting Information S1), the systems with L = 22.5 m had smaller percolation thresholds compared to those with L = 30, 50, and 70 m (Figure S4 in Supporting Information S1 To see the fits to the rest of the cases, see FiguresS2 and S3in Supporting Information S1.