A simplified relationship between the zero-percolation threshold and fracture set properties

Percolation analysis is an efficient way of evaluating the connectivity of discrete fracture networks. Except for very simple cases, it is not feasible to use analytical approaches to find the percolation threshold of a discrete fracture network. The most commonly used percolation threshold corresponds to the occurrence of percolation on average for the set of parameters (p50), which is not adequate for applications in which a high confidence in the percolation threshold is required. This study investigates the direct relationships between the percolation threshold at low probability (p0, referred to as zero-percolation threshold) and the properties of fracture networks with one set of fractures (fractures with similar orientations) in two-dimensional domains. A generalized non-linear multivariate relationship between p0 and fracture network parameters is established based on connectivity assessments of a significant number of numerical simulations of fracture networks. A feature of this relationship is the invariant shape of marginal relationships. A comparison study with an analytical solution and applications in both synthetic and real fracture networks shows that the derived relationship performs well in fracture networks of different sizes and orientations. A significant benefit of this relationship is that, when an analytical solution is not available, it can provide fast and reliable connectivity statistics of fracture networks based only on fracture parameters.

1 Introduction

Discrete fracture networks (DFN) are widely used to model fracture systems in rock masses and reservoirs for flow analysis (Dogan, 2023; Kolyukhin, 2022; Liu et al., 2019). In three-dimensional DFNs, fractures are represented using simplified geometric shapes such as polygons, rectangles, disks, and others. In two-dimensional DFNs, fractures are simplified as line segments (Dong et al., 2020). DFN models are generally considered to be more realistic representations of the fracture systems and they are more amenable to integrating geological data into the flow model (Dong et al., 2018a). Fracture connectivity analysis is a significant component in the DFN approach to assessing flow behaviour (Alghalandis et al., 2015; Einstein and Locsin, 2012) as the conductivities of fractures are generally orders of magnitude greater than those of the surrounding porous matrix (Thovert et al., 2017). A good understanding of the connectivity of fracture networks is essential for many applications such as oil and gas recovery, geothermal energy exploitation, hydrology and groundwater engineering, and geological storage of radioactive wastes.

Percolation theory (Jafari and Babadagli, 2013; Khamforoush and Shams, 2007; Or et al., 2023; Sun et al., 2023) provides a basis for describing and quantifying the connectivity of geometrically complex systems (Xu et al., 2007) such as fracture networks and this is reflected in many studies reported in the literature that use percolation theory in the connectivity analysis of fracture systems (Barker, 2018; de Dreuzy et al., 2000; Dong et al., 2022; Khamforoush and Shams, 2007; Manzocchi, 2002; Masihi and King, 2007; Mourzenko et al., 2012). Percolation describes the phenomenon in which there is at least one domain-spanning pathway in a physical system (Tang et al., 2022; Yao et al., 2020; Yi and Tawerghi, 2009). Percolation in a DFN involves at least one cluster of connected fractures that spans the reservoir (McKenna et al., 2020) or rock mass (one cluster refers to a series of fractures that intersect each other). In this context, one of the most important characteristics of a fracture network is whether or not it percolates (Bour and Davy, 1998; Bour and Davy, 1997) and the percolation threshold is commonly used to quantify the critical value of connectivity at which the network percolates (Khamforoush and Shams, 2007; Manzocchi et al., 2023; Walsh and Manzocchi, 2021). Using this definition, the permeability of a fracture network is zero if the connectivity value is less than the percolation threshold (Mourzenko et al., 2005).

The percolation threshold of a DFN is a property that depends on the parameters of the fracture system (Mourzenko et al., 2005). Features previously used to characterise percolation in DFNs include the dimensionless density derived from the excluded volume (Barker, 2018; de Dreuzy et al., 2000; Khamforoush et al., 2008; Mourzenko et al., 2012), fractal dimensions (Jafari and Babadagli, 2013; Jafari and Babadagli, 2009; Zhao et al., 2016), topological connectivity measures (Manzocchi, 2002), fracture clustering (Manzocchi, 2002), and the average number of intersections per fracture (Manzocchi, 2002). These indirect characteristics of DFN models are derived from direct fracture network parameters, such as the number of fractures, fracture locations, sizes, and orientations, which have a joint effect on the occurrence of a percolating network (Jafari and Babadagli, 2009). For example, if the fracture size is kept constant, an increase in the number of fractures will result in a higher fracture density, which in turn will increase the probability of a connected domain (Shokri et al., 2016).

There are many published studies on the percolation of DFN models, with different focuses on different aspects of the problem. For fracture locations, DFN models in these studies cover both the Poisson (homogeneous) distribution (Barker, 2018; Bour and Davy, 1997; de Dreuzy et al., 2000; Huseby and Thovert, 1997; Jafari and Babadagli, 2013; Khamforoush and Shams, 2007; Mourzenko et al., 2005; Robinson, 1983; Thovert et al., 2017; Zhao et al., 2009; Zhao et al., 2016) and non-homogeneous (i.e., spatially correlated) distributions (Manzocchi, 2002; Mourzenko et al., 2012). For fracture sizes, some DFNs use the monodisperse model, which means that the shape and size of every fracture are identical (Jafari and Babadagli, 2013; Khamforoush and Shams, 2007; Khamforoush et al., 2008; Manzocchi, 2002; Mourzenko et al., 2012; Robinson, 1983). This makes the percolation study relatively simple using numerical simulation. Others use the polydisperse model in which the sizes (Bour and Davy, 1997; Huseby and Thovert, 1997; Mourzenko et al., 2004; Mourzenko et al., 2005; Thovert et al., 2017; Zhao et al., 2016) and shapes (Barker, 2018; Thovert et al., 2017) of fractures are different. Commonly used fracture size distributions include the power law function (Mourzenko et al., 2004; Mourzenko et al., 2005; Zhao et al., 2016), exponential distribution (Catapano et al., 2023; Dowd et al., 2007; Fadakar Alghalandis, 2017; Xu et al., 2007; Zhu et al., 2022) and uniform distribution (Huseby and Thovert, 1997). For fracture orientations, many DFNs use isotropic models (uniform and random) (Barker, 2018; Bour and Davy, 1997; Charlaix et al., 1984; de Dreuzy et al., 2000; Huseby and Thovert, 1997; Jafari and Babadagli, 2013; Khamforoush and Shams, 2007; Khamforoush et al., 2008; Mourzenko et al., 2004; Mourzenko et al., 2012; Mourzenko et al., 2005; Robinson, 1983; Thovert et al., 2017; Yi and Tawerghi, 2009; Zhao et al., 2016) but anisotropic (with a single preferential orientation or several preferential orientations) models are also commonly used (Balberg and Binenbaum, 1983; Khamforoush and Shams, 2007; Khamforoush et al., 2008; Manzocchi, 2002). The Fisher distribution (Khamforoush and Shams, 2007; Xu and Dowd, 2010) is the most commonly used type of distribution for three-dimensional fracture networks, while the von Mises distribution (Xu and Dowd, 2010) is the most commonly used for two-dimensional fracture networks. Physically, fractures related to tectonic movements are, in general, anisotropic (e.g., conjugate fractures generated around the maximum principal compressive stress (Zhao and Hou, 2017), while fractures associated with other causes, such as diagenesis, are typically isotropic (Dong et al., 2018b).

A common method used to obtain the percolation threshold of a DFN is first to calculate some indirect characteristic parameters of the fracture network and then evaluate the percolation threshold on the basis of these parameters. However, DFNs with the same indirect characteristic parameters may have quite different direct geometrical parameters (e.g., number of fractures, fracture size, and orientation). Unfortunately, these geometrical parameters dictate the fracture connectivity and hence the percolation threshold (Dong et al., 2019). For example, the two fracture networks in Fig. 1 have the same number of fractures, identical fracture lengths, and box-counting fractal dimensions, but the network in Fig. 1a percolates between sides A and B, while the other (Fig. 1b) does not percolate. The different orientations of these two fracture models result in different percolation characteristics. In this case, the box-counting fractal dimension provides a good measure of the complexity of the system but it ignores the effect of the preferential orientation of a fracture network. Although the indirect approach can simplify the evaluation of the percolation threshold of a fracture network, it may sometimes produce misleading results. In addition, most percolation thresholds based on the excluded volume method correspond to the occurrence of percolation on average (Barker, 2018; Yi and Tawerghi, 2009) (i.e., with 50 % probability, p50) due to the stochastic nature of fracture networks (cf. Sect. 2.2). The level of confidence in such thresholds may not be sufficient for some applications. For example, for underground radioactive waste storage facilities, selecting storage sites that minimize the potential number of connected pathways to the biosphere is important. In this case, a low probability percolation threshold (p0, cf., Sect. 2.2) named as zero-percolation threshold is more important for the connectivity analysis of the fracture systems (Dong et al., 2019).

Figure 1 Schematic diagram of two DFN models with different percolation features. (a) Percolated DFN; (b) Un-percolated DFN.

In practical engineering applications, it would be less problematic to use the direct relationship between the percolation threshold and fracture network parameters for connectivity assessments instead of resolving the problem by numerical simulations in every case. In general, it is not possible to establish such a relationship analytically due to the complexity of DFNs. In this work, Monte Carlo simulations of the number of fractures, n, fracture size (length), L ,and fracture orientation, ϕ are used to establish this relationship for two-dimensional DFN models, the fractures are represented by vectors representing line segments (Dong et al., 2020). In particular, fracture locations follow the Poisson distribution, fracture lengths follow the exponential distribution f(L|λ) (λ: the exponential distribution parameter whose reciprocal 1/λ represents the mean fracture length $\stackrel{\mathrm{‾}}{L}$) and fracture orientations follow the von-Mises distribution f(ϕ|μ,κ), where λ,μ and κ are their corresponding distribution parameters. μ represents the mean/main fracture orientation, while κ (concentration parameter) quantifies the directional clustering degree of fractures in the von Mises distribution (see Sect. 2.1). The zero-percolation threshold (p0) equation $L_{\mathrm{t}}=f\left(n,\mathit{\mu },\mathit{\kappa }\right)$ was established by analysing results from an extensive set of numerical simulations for fracture networks with one set of fractures (fractures that exhibit similar orientations) (Ali and Jakobsen, 2011; Zeng et al., 2022) (see Sect. 3.1–3.2). Besides DFN with exponential fracture length, given the widespread adoption of lognormal distribution f(L|μσ) in characterizing fracture length distributions within DFN, it is imperative to explore the implications of this distribution on the phenomenon of zero-percolation. Here, μ and σ are parameters of the lognormal distribution, where μ is the mean of the natural logarithm of fracture lengths and σ is the standard deviation of the natural logarithm of fracture lengths (derived from the mean and standard deviation of fracture lengths). Consequently, this paper extends its investigation beyond DFNs employing exponential distribution for fracture lengths to encompass those utilizing lognormal distribution, as detailed in Sect. 3.3–3.4.

The relationship was established by using a non-linear, multivariate fitting method for a relationship with invariant shapes of marginal functions; this is demonstrated by a simple example in Sect. 2.2 and 2.3. The verification of the derived equations for zero-percolation underwent a comprehensive series of tests in Sect. 4.

2 Principle of the mathematical method

To mitigate ambiguity, the mathematical approach will utilise DFN with exponential distributions as an instance to elucidate the fundamental principles. Fracture parameters of a DFN model are described in Sect. 2.1. The percolation and percolation threshold of a DFN are described in Sect. 2.2. In Sect. 2.3, the non-linear fitting method is illustrated using a simple percolation example.

2.1 Discrete fracture networks

DFN modelling is a stochastic simulation method that uses marked point processes (MPP) (Dong et al., 2018b) in which fracture location (x,y) is modelled by a point process (Fig. 2a) following a Poisson, non-homogeneous cluster or Cox point process (Mardia et al., 2007); and fracture properties (such as length L, orientation ϕ) are modelled at each point by marks (Fig. 2b) following their respective probability distribution functions (e.g., f(L) and f(ϕ)) (Dong et al., 2018b; Fadakar Alghalandis, 2017; Xu and Dowd, 2010). To simulate a set of n fractures with similar orientations, the location of a fracture (Fig. 2a) is generated first followed by the generation of the associated marks (Fig. 2b). Subsequently, these procedures are iterated n times to culminate in the ultimate implementation of the DFN (Fig. 2c).

Figure 2 Schematic diagram of a two-dimensional DFN realization. (a) Randomly generated fracture location; (b) Fracture properties (length, orientation) are generated from their probability distributions; (c) Repeat process (a) and (b) to generate the entire fracture network to account for the number of fractures in each fracture set as well as the number of fracture sets. (a) Point process; (b) Mark generation; (c) Discrete fracture network.

For 2D DFN models, fracture network parameters include the number of fracture sets, number of fractures n in each fracture set, fracture size distribution f(L) and fracture orientation distribution f(ϕ). The study area used in this work is 100 m × 100 m so $P\mathrm{20}=n_i\times {\mathrm{10}}^{-\mathrm{4}}\left({\mathrm{m}}^{-\mathrm{2}}\right)$. P20 (2D fracture density) represents the number of fractures per unit area. Here, P20 is the fracture number per 2D unit area (Khamforoush et al., 2008). For fracture length, a fixed size L can be used, or the following exponential distribution (Xu and Dowd, 2010) is commonly used:

$$\begin{array}{}\text{(1)}& f\left(L\mathrm{|}\mathit{\lambda }\right)=\left\{\begin{array}{ll}\mathit{\lambda }{e}^{-\mathit{\lambda }L}& L\ge \mathrm{0}\\ \mathrm{0}& L<\mathrm{0}\end{array}\right.,\end{array}$$

where λ is the distribution parameter and therefore the average length, which in this case is $\stackrel{\mathrm{‾}}{L}=\mathrm{1}/\mathit{\lambda }$. For fracture orientation, a von-Mises distribution (Fadakar Alghalandis, 2017; Xu and Dowd, 2010) is commonly used for 2D applications:

$$\begin{array}{}\text{(2)}& f\left(\mathit{\varphi }\mathrm{|}\mathit{\mu },\mathit{\kappa }\right)={\displaystyle \frac{e^{\mathit{\kappa }\cos (\mathit{\varphi }-\mathit{\mu })}}{\mathrm{2}\mathit{\pi }{I}_{\mathrm{0}}\left(\mathit{\kappa }\right)}},\end{array}$$

where μ and κ are the distribution parameters, I₀ is the modified Bessel function defined as: $I_{\mathrm{0}}\left(\mathit{\kappa }\right)={\sum }_{i=\mathrm{0}}^{+\mathrm{\infty }}\frac{{\mathit{\kappa }}^{\mathrm{2}i}}{{\mathrm{2}}^{\mathrm{2}i}(i\mathrm{!}{)}^{\mathrm{2}}}$.

μ and $\mathrm{1}/\mathit{\kappa }$ are analogous to the mean and variance of a Gaussian distribution. μ represents the main fracture orientation. For example, in both Fig. 3a and b, the main orientation is NE-SW (Fig. 3e and f) so $\mathit{\mu }=\mathit{\pi }/\mathrm{4}$, while in Fig. 3c and d, $\mathit{\mu }=\mathit{\pi }/\mathrm{2}$ (Fig. 3g and h), following the common practice in geotechnical applications of measuring the bearing angle from the North. κ is a measure of concentration (reciprocal of the measure of dispersion). A comparison of Fig. 3a and c with b and d shows that the dispersion of the fracture orientation increases as $\mathrm{1}/\mathit{\kappa }$ increases. When κ=0, the distribution of fracture orientations is completely random. The DFN models were generated by a Matlab code based on the open-source toolbox ADFNE (Fadakar Alghalandis, 2017).

Figure 3 DFN models showing different fracture orientations with different μ, κ, together with their corresponding rose diagrams. Fractures of the same colour in the DFN model are in the same cluster, while fractures in black are isolated ones. (a) μ=π/4, κ= 4; (b) μ=π/4, κ= 24; (c) μ=π/2, κ= 4; (d) μ=π/2, κ= 24; (e) μ=π/4, κ= 4; (f) μ=π/4, κ= 24; (g) μ=π/2, κ= 4; (h) μ=π/2, κ= 24.

2.2 Percolation threshold of DFN models

Percolation in a DFN means there is at least one cluster of fractures spanning the system (rock mass or reservoir) that allows the fluid to permeate from one side to the other, as shown in Fig. 4b, while Fig. 4a shows a non-percolating DFN. Obviously, one can easily conclude that an increase in fracture number or fracture size can lead to the percolation of the system at some stage. In reality, the probability of percolation is a function of many different factors related to the DFN parameters (Khamforoush et al., 2008), of which fracture density (such as P20 and P21 for 2D and P30 and P32 for 3D applications) is the most critical. Here, P21 is the fracture length per 2D unit area, while P30 and P32 are the fracture number and area per 3D unit volume, respectively.

Figure 4 Schematic diagram of percolation in DFN models. (a) Un-percolated DFN; (b) Percolated DFN.

For simple fracture networks, the percolation thresholds can be found analytically. However, for most fracture networks, approaches such as Monte Carlo (MC) simulation are required (Yi and Tawerghi, 2009). Due to the random nature of a DFN model, its corresponding percolation status is also stochastic in nature. If N independent simulations of a DFN are repeated for a group of parameters $\left(n,L,\mathit{\varphi }\right)$, resulting in N_p number of cases where the fracture network percolates, then the percolation probability corresponding to the parameter set is $P=N_p/N$ (Barker, 2018; Yi and Tawerghi, 2009). In this study, p0, p50, and p100 represent the percolation thresholds of the DFN corresponding to 0 %, 50 %, and 100 % percolation, respectively. p0 (zero-percolation threshold) is defined as the critical value of a fracture parameter below which network connectivity is completely lost, ensuring absolute impermeability. In a fracture network with two variables, fracture number (n) and fracture length (L), p0 may be represented either by a sufficiently small n (regardless of L) or by a minimum L when n is fixed.

As illustrated in Fig. 5, percolation probability typically transitions between P=0 (no percolation) and P=1 (full percolation) through an intermediate band. The choice of percolation threshold depends on the application. For example, in subterranean radioactive waste repositories, where preventing any potential connectivity to aquifers is critical (Wei et al., 2017; Yi and Tawerghi, 2009), p0 is adopted to ensure absolute impermeability. Conversely, for applications such as unconventional gas extraction or enhanced geothermal systems, near-complete connectivity is desired, and the threshold may correspond to P≈1. p50 is often adopted in resource extraction (e.g., hydrocarbon reservoirs) to balance economic viability and manageable risk. Noted in stochastic systems, P= 0 % and 100 % may not be strictly possible and therefore the definition used here means the probability calculated by $N_{\mathrm{p}}/N_T$ using a reasonable number N_T( = 20 in this study). For stochastic systems, Fig. 5 shows a worked example, in which the fracture orientations follow a von-Mises distribution ($\mathit{\mu }=\mathrm{90}\mathit{°}\mathit{\kappa }=\mathrm{24}$) and fracture lengths are identical for each DFN. Twenty MC simulations (N_T=20) were conducted using pairs of parameters ($n=\mathrm{20},\mathrm{30},\mathrm{40},\mathrm{\dots },\mathrm{250}$ and $L=\mathrm{0.06},\mathrm{0.08},\mathrm{0.1},\mathrm{\dots }\mathrm{0.8}$). The percolation probability, P, calculated from the simulations, is shown in Fig. 5, where the horizontal axes correspond to the number of fractures n and the fracture length L, respectively, and the z axis is the percolation probability of the corresponding DFN model.

Figure 5 Percolation probability P versus fracture number n and fracture length L for DFN models with random orientations and identical lengths. Each P value is obtained from 20 Monte Carlo simulations for a given (n, L) pair. For each n, the fracture length threshold corresponding to the p0? zero-percolation point is determined, and the average over 20 repeats is taken as the statistical p0 threshold.

Figure 6a presents data points in the vicinity of the p0 percolation threshold. For each fracture number (n), the statistical mean of the fracture length threshold (L_t) corresponding to p0 percolation is determined from 20 independent Monte Carlo (MC) simulations, following the same procedure as in Fig. 5 with the same parameters ($\mathit{\mu }=\mathrm{90}\mathit{°},\mathit{\kappa }=\mathrm{24}$). The mean values are depicted as circles, with half error bars indicating the mean ± two standard deviations in Fig. 6a. To further approach the true threshold values, an additional 20 repetitions of the Fig. 6a procedure are performed, providing a statistically robust estimate that closely approximates the true value. In order to reduce the computation cost, there are a few differences in fracture number n and length L compared with those in Fig. 5, i.e., $n=\mathrm{round}\left({\mathrm{10}}^j\right)=\mathrm{22},\mathrm{32},\mathrm{45},\mathrm{\dots },\mathrm{251},j=\mathrm{1.35},\mathrm{1.5},\mathrm{1.65},\mathrm{\dots },\mathrm{2.4},L={\mathrm{10}}^k=\mathrm{0.0316},\mathrm{0.0322},\mathrm{0.0327},\mathrm{\dots },\mathrm{1},k=-\mathrm{1.5},-\mathrm{1.4925},-\mathrm{1.485},\mathrm{\dots },\mathrm{0}$. The corresponding non-linear percolation threshold curve based on least-squares regression is $L_{\mathrm{t}}=\mathrm{1.7518}\times n^{-\mathrm{0.4308}}$ with a correlation coefficient of 0.9288. When the specific number n is considered, an increase in L results in percolation. In this context, we hold the fracture number n to determine the fracture length threshold (denoted as L_t) for zero probability percolation, hence the utilisation of L_t and n. Conversely, if the fracture length is fixed, Land n_t will be employed. This curve defines the percolation threshold in terms of parameter pairs of (n,L_t). The DFN corresponding to any combination of parameters below this curve, namely L<L_t, will have a percolation probability of 0. To assess the uncertainty of this relationship, 20 groups of MC simulations as used in Fig. 6a were employed to obtain the non-linear relationship with statistical dispersion and the results are shown in Fig. 6b. The circles correspond to the averages of 20 groups of 20 MC simulations (Fig. 6a) and the error bars represent two times the corresponding standard deviations. The fitted non-linear relationship based on least-squares regression of the average values is $L_{\mathrm{t}}=\mathrm{1.592}\times n^{-\mathrm{0.403}}$ with a correlation coefficient of 0.993.

Figure 6 Percolation threshold of (n, L_t) for the example DFN model. (a) L_t vs. n from 20 MC simulations in Fig. 5, with circles showing means and half error bars indicating twice the standard deviation. (b) (n, L_t) corresponding to different Δ, where each point is the average of 20 repetitions in (a).

Fracture orientations will also affect the percolation threshold curves. As the example to demonstrate these effects, the von-Mises distribution is used to describe the distribution of fracture orientations for the DFN model used above. In this case, the fracture network parameters are ($n,L,\mathit{\mu },\mathit{\kappa }$). To simplify the demonstration, the concentration parameter κ is set to 24, similar to those shown in Fig. 3b and d;  μ is set to values from 90 to 0° in 5° decrements. As the horizontal percolation is of interest here, to simplify the comparison, μ is transformed to an angle measured from the horizontal direction, i.e., $\mathrm{\Delta }=\left|\mathit{\mu }-\mathrm{90}\right|$(Relative orientation angle, indicating the deviation angle between the dominant fracture direction μ and the horizontal direction, used for simplified horizontal percolation analysis), hence, $\mathrm{\Delta }=\mathrm{0}\mathit{°},\mathrm{5}\mathit{°},\mathrm{\dots },\mathrm{90}\mathit{°}$. The above curve fitting process was repeated and some results are shown in Fig. 6b. As Δ decreases, the percolation threshold decreases. This is consistent with the fact that lower Δ will increase the connectivity in the horizontal direction between the left and the right sides; percolation therefore requires a shorter fracture length, and so the percolation threshold curve decreases.

2.3 Non-linear relationship for the percolation threshold

The general shapes of the percolation curves in Fig. 6c can be considered to have similar shapes, which are represented as coloured curves in Fig. 7. To understand further the relationship between L_t, n and Δ, average values of 20 MC simulations are shown in Fig. 7. Coloured curves refer to slices of different n and Δ. The left colour curves are the variation of L_t and n with different values of Δ, and the right coloured curves are the variations of L_t and Δ with different values of n. Clearly, the threshold fracture length increases as the number of fractures decreases and the relative orientation Δ increases. This is because higher fracture density leads to greater percolation probability, whereas a lower number of fractures requires longer fractures to maintain the same percolation probability.

Figure 7 Simulated percolation threshold L_t vs. (n,Δ). For each pair of (n,Δ), 20 MC simulations are implemented to obtained the average L_t. Coloured curves are lines corresponding to slices of different n and Δ.

From the results in Fig. 7, $L_{\mathrm{t}}=f\left(n,\mathrm{\Delta }\right)$ is non-linear. To establish this relationship, the variation of L_t with n for different values of Δ is examined first, i.e., L_t=f₁(n)|_Δ, followed by assessing the influence of Δ on the derived f₁(n) relationship. Note that at the second stage, cos Δ is used instead of Δ as it is more relevant to the quantification of the fracture projection length in the horizontal direction ($\mathrm{len}/\cos \mathrm{\Delta })$.

Based on the simulation results discussed above, Eq. (3) is considered an appropriate fit to f₁(n):

$$\begin{array}{}\text{(3)}& L_{\mathrm{t}}=f_{\mathrm{1}}\left(n\right)=a{n}^b,\end{array}$$

where a and b are parameters to be determined in the fitting process. The correlation coefficients for all curves for different Δ values (0, 5 … 85, 90°) are 0.9985, 0.9895, … , 0.9949, respectively. The high correlation coefficients (> 0.96) confirm the suitability of using Eq. (3) to represent f₁(n).

The relationships between a,b and cos Δ, shown in Fig. 8, suggest linear relationships. Least squares regression was used to obtain the following equations:

$$\begin{array}{}\text{(4)}& {\displaystyle }a=a_{\mathrm{1}}{\cos }^{\mathrm{2}}\mathrm{\Delta }+a_{\mathrm{2}}\cos \mathrm{\Delta }+a_{\mathrm{3}},\text{(5)}& {\displaystyle }b=b_{\mathrm{1}}{\cos }^{\mathrm{2}}\mathrm{\Delta }+b_{\mathrm{2}}\cos \mathrm{\Delta }+b_{\mathrm{3}},\end{array}$$

with correlation coefficients of 0.9919 and 0.9712, respectively; a₁, a₂, a₃, b₁, b₂, b₃ are 2.726, −5.304, 2.887, −0.2552, 0.6134, and −0.5724, respectively.

Figure 8 Influence of Δ on a and b. (a) a vs. cos Δ; (b) b vs. cos Δ.

By incorporating Eqs. (4) and (5) into Eq. (3) the final form of the expression of L_t in terms of the fracture network parameters is obtained, as shown in Eq. (6). This form is then used directly in a bivariate least squares fitting using the Levenberg-Marquardt algorithm (Ngia and Sjoberg, 2000), an optimal search technique for multivariate non-linear curve fitting. The original values of parameters shown in Eqs. (4) and (5) are used as initial inputs to the optimisation algorithm to improve computational efficiency and accuracy and the final derived parameters in this case are (a₁, a₂, a₃, b₁, b₂, b₃) = (2.4757, −4.9064, 2.7359, −0.1841, 0.5097, −0.5336). This set of values should be a more accurate reflection of the bivariate relationship than the values obtained in the two separate consecutive steps described above. The ${\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}$ parameter (mean fracture length threshold) quantifies the orientation-dependent critical length for percolation, defined mathematically as:

$$\begin{array}{}\text{(6)}& {\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=\left(a_{\mathrm{1}}{\cos }^{\mathrm{2}}\mathrm{\Delta }+a_{\mathrm{2}}\cos \mathrm{\Delta }+a_{\mathrm{3}}\right)n^{\left(b_{\mathrm{1}}{\cos }^{\mathrm{2}}\mathrm{\Delta }+b_{\mathrm{2}}\cos \mathrm{\Delta }+b_{\mathrm{3}}\right)},\end{array}$$

The final fitted surface is shown in Fig. 9a. The points are the average values of 20 groups of MC simulation results shown in Fig. 7. The suitability of the chosen functional form (Eq. 6) is confirmed by the fact that almost all the points are on the fitted surface. The plot in Fig. 9b of simulated values of L_t against those predicted by Eq. (6) gives an extremely high correlation coefficient of nearly 1. It is also encouraging that, on visual inspection, the fitted curve is conditionally unbiased.

Figure 9 Bivariate percolation threshold fitting. (a) Fitted percolation threshold surface and simulation data; (b) Cross plot between simulated values of L_t and those predicted by Eq. (6).

Although this workflow is useful for multivariate non-linear fitting problems in which marginal relationships are of invariant shape, it should be noted that difficulties may arise for very high dimensions (Dong et al., 2016).

The process described above can be summarised as an approach for fitting multiple variables. This approach starts by fitting a hypothetical relationship between L_t and n is initially fitted. Then, a new variable Δ is added by analysing relationships with the parameters in the hypothetical relationship model. The parameters in the hypothetical relationship model are then replaced by expressions of the newly added variable. Ultimately, the relationship between L_t and (n, Δ) can be obtained. This approach will be applied in the relationship fitting of Sect. 3.2, where the independent variables are (ncos Δκ).

3 Percolation analysis of DFN models

3.1 Experiment design for DFN with exponential fracture lengths

In the example used above, the lengths and orientations were identical for all fractures in a fracture network, which was not generally the case in practical applications. The relationship described above could be made more useful by extending it to cover realistic fracture networks. The following numerical experiments were all implemented on a dimensionless unit square (1×1). Based on previously published work (Dong et al., 2018c; Xu et al., 2007), the lengths of rock fractures could generally be modelled by an exponential or lognormal distribution. In this work, the exponential distribution was used, and therefore the average length $\stackrel{\mathrm{‾}}{L}$ is equal to $\mathrm{1}/\mathit{\lambda }$ where λ was the distribution parameter. The von-Mises distribution (Eq. 2) was used for fracture orientation.

There are now three independent variables $(n,\mathrm{\Delta },\mathit{\kappa })$ and the aim is to establish the relationship ${\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=f(n,\mathrm{\Delta },\mathit{\kappa })$. To simulate percolation states similar to those shown in Fig. 5, DFNs corresponding to a combination of $\mathrm{8}\times \mathrm{18}\times \mathrm{13}\times \mathrm{200}$ (374 400) sets of variables were simulated and analysed, with each case simulated independently 20 times. The number of changes explored for each variable are listed in Table 1.

Table 1 Parameters of DFNs in Sect. 3.2.

For each pair of (Δ,κ), 20 independent realisations of DFNs with different n and $\stackrel{\mathrm{‾}}{L}$ were generated to obtain the percolation threshold curves ${\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=f_{\mathrm{1}}\left(n\right)$. These 20 MC simulations are used to calculate the percolation probability to obtain the points $\left(n,{\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}\right)$, as shown in Fig. 10. The points are the average values of 20 groups of 20 MC simulations and the error bars represent two times of the corresponding standard deviation. The relationship in Eq. (3) was used again for ${\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=f_{\mathrm{1}}\left(n\right){\mathrm{|}}_{\mathrm{\Delta },\mathit{\kappa }}$. A comparison of Figs. 6b and 10 indicates that the standard deviations in this case are much larger, which is expected due to the variability in the lengths of fractures generated in simulations. Note that the uncertainty (reflected by the size of the error bar) increases as the number of fractures, n, decreases.

Figure 10 Variation of the percolation threshold as a function of n for κ=3 and Δ=0.

3.2 Determining percolation threshold equation for DFN with exponential fracture length

Equation (6) is only for a specific parameter κ. Establishing the full relationship, ${\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=f(n,\cos \mathrm{\Delta },\mathit{\kappa })$, requires the relationships between (a₁, a₂, a₃, b₁, b₂, b₃) in Eq. (6) and κ for different DFNs. The regression results for these parameters for different κ values are shown in Table 2, which provides the input data for the final non-linear fitting of the percolation threshold function. The correlation coefficients for each of the regressions in the table are all greater than 0.98, which ensures the suitability of the derived relationships.

Table 2 Regression parameters (a₁, a₂, a₃, b₁, b₂, b₃) for different values of κ.

Figure 11 Relationship between κ and fitted parameters (a₁, a₂, a₃, b₁, b₂, b₃). (a) κ vs. a₁; (b) κ vs. a₂; (c) κ vs. a₃; (d) κ vs. b₁; (e) κ vs. b₂; (f) κ vs. b₃.

The relationships between a₁ and κ, between a₂ and κ and between a₃ and κ are linear as described by Eqs. (7)–(9); the relationship between b₁ and κ, between b₂ and κ and between b₃ and κ are quadratic as shown in Eqs. (10)–(12). The correlation coefficients of this set of regression curves are all greater than 0.98. Overall, these variables display a clear and strong relationship that can be described by an appropriate functional form. Table 3 lists the constants in Eqs. (7)–(12) obtained by least squares regression.

$$\begin{array}{}\text{(7)}& {\displaystyle }{a}_{\mathrm{1}}=a_{\mathrm{11}}\mathit{\kappa }+a_{\mathrm{12}},\text{(8)}& {\displaystyle }{a}_{\mathrm{2}}=a_{\mathrm{21}}\mathit{\kappa }+a_{\mathrm{22}},\text{(9)}& {\displaystyle }{a}_{\mathrm{3}}=a_{\mathrm{31}}\mathit{\kappa }+a_{\mathrm{32}},\text{(10)}& {\displaystyle }{b}_{\mathrm{1}}=b_{\mathrm{11}}\mathit{\kappa }+b_{\mathrm{12}}{\mathit{\kappa }}^{\mathrm{2}}+b_{\mathrm{13}},\text{(11)}& {\displaystyle }{b}_{\mathrm{2}}=b_{\mathrm{21}}\mathit{\kappa }+b_{\mathrm{22}}{\mathit{\kappa }}^{\mathrm{2}}+b_{\mathrm{23}},\text{(12)}& {\displaystyle }{b}_{\mathrm{3}}=b_{\mathrm{31}}\mathit{\kappa }+b_{\mathrm{32}}{\mathit{\kappa }}^{\mathrm{2}}+b_{\mathrm{33}},\end{array}$$

where a₁₁, a₁₂, a₂₁, a₂₂, a₃₁, a₃₂, b₁₁, b₁₂, b₁₃, b₂₁, b₂₂, b₂₃, b₃₁, b₃₂, b₃₃ are parameters.

Table 3 Parameters of Eqs. (7)–(12).

Finally, the combined percolation equation ${\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=f(n,\cos \mathrm{\Delta },\mathit{\kappa })$ can be obtained, as shown in Eq. (13). The correlation coefficient is again nearly 1 (0.99).

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}& =\left(\right(a_{\mathrm{11}}\mathit{\kappa }+a_{\mathrm{12}})x^{\mathrm{2}}+(a_{\mathrm{21}}\mathit{\kappa }+a_{\mathrm{22}})x\\ & +(a_{\mathrm{31}}\mathit{\kappa }+a_{\mathrm{32}}))n^{\left(\begin{array}{l}(b_{\mathrm{11}}\mathit{\kappa }+b_{\mathrm{12}}{\mathit{\kappa }}^{\mathrm{2}}+b_{\mathrm{13}})x^{\mathrm{2}}\\ +(b_{\mathrm{21}}\mathit{\kappa }+b_{\mathrm{22}}{\mathit{\kappa }}^{\mathrm{2}}+b_{\mathrm{23}})x+b_{\mathrm{31}}\mathit{\kappa }\\ +b_{\mathrm{32}}{\mathit{\kappa }}^{\mathrm{2}}+b_{\mathrm{33}}\end{array}\right)},\end{array}\end{array}$$

where x=cos Δ. x is a derived variable that optimizes the equation structure.

Again, these parameters are used as the inputs for the final multivariate least squares optimisation based on the Levenberg-Marquardt algorithm using all the simulation results. The final optimised values of the parameters in Eq. (13) are shown in Table 4. These values are similar to the initial parameter values obtained by the step-wise fitting process described above but they have been refined by global optimisation. The correlation coefficient between the prediction and simulation values based on the initial parameters (Table 3) is only 0.43 due to the error propagation in the step-wise fitting process. After global optimisation, the correlation coefficient increases significantly to nearly 1 (0.99) based on the values listed in Table 4.

Table 4 Parameters of the percolation equation in terms of fracture properties.

To visualize the relationships in Eq. (13), several surfaces of ${\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}$ vs. (n,Δ) corresponding to different values κ (4, 18, 42 and 75) are shown in Fig. 12. In general, higher κ values correspond to higher percolation threshold values. This is because higher κ values correspond to lower variation of fracture orientations, which leads to lower probabilities of fracture intersections. Consequently, this reduces the connectivity of the fracture network and hence longer fractures are needed to reach percolation.

Figure 12 Extracted surfaces from the final percolation threshold equation. (a) κ=4; (b) κ=18; (c) κ=42; (d) κ=75.

Equation (13) is derived for the region of a dimensionless unit square, the result is expected to be applicable to areas at different scales (y×y). For these cases, the scaled percolation threshold ${\stackrel{\mathrm{‾}}{L}}_{ts}$ will be used (Eq. 14) revised from Eq. (13). If the average fracture length of a fracture network $\stackrel{\mathrm{‾}}{L}\ge {\stackrel{\mathrm{‾}}{L}}_{\mathrm{ts}}$, there is more probability reach percolation.

$$\begin{array}{}\text{(14)}& {\stackrel{\mathrm{‾}}{L}}_{\mathrm{ts}}={\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}\times y,\end{array}$$

3.3 Design of experiments for DFN with lognormal fracture lengths

Unlike the previous example, the fracture lengths in this section are modeled with a lognormal distribution. The mean fracture length $\stackrel{\mathrm{‾}}{L}$ and standard deviation ν, which quantifies the dispersion of fracture lengths under the lognormal distribution, are used to derive the distribution parameters μ and σ. These parameters, along with the probability density function, are calculated according to Eqs. (15)–(17). Five independent variables $(\stackrel{\mathrm{‾}}{L},\mathit{\nu },\mathrm{\Delta }\mathit{\kappa })$ are considered with the aim of establishing the relationship $n_{\mathrm{t}}=f\left(\stackrel{\mathrm{‾}}{L},\mathit{\nu },\mathrm{\Delta },\mathit{\kappa }\right)$. In this context, n_t represents the fracture number threshold at which the percolation threshold may be reached at a low probability (p0) in DFNs characterized by the parameters $(\stackrel{\mathrm{‾}}{L},\mathit{\nu },\mathrm{\Delta }\mathit{\kappa })$. In Sect. 3.1, the exponential distribution of fracture length is defined by a single parameter. Consequently, the fracture length is selected to determine the threshold corresponding to p0. Given that the fracture length is governed by two parameters $(\stackrel{\mathrm{‾}}{L},\mathit{\nu })$, the parameters corresponding to the fracture length distribution are not selected. Instead, the fracture number is chosen.

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathit{\mu }=\log ({\stackrel{\mathrm{‾}}{L}}^{\mathrm{2}}/\sqrt{\mathit{\nu }+{\stackrel{\mathrm{‾}}{L}}^{\mathrm{2}}}),\text{(16)}& {\displaystyle }\mathit{\sigma }=\sqrt{\log (\mathit{\nu }/{\stackrel{\mathrm{‾}}{L}}^{\mathrm{2}}+\mathrm{1})},\text{(17)}& {\displaystyle }f\left(\stackrel{\mathrm{‾}}{L}\right|\mathit{\mu },\mathit{\sigma })={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{L}\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}}{\mathrm{e}}^{\frac{-(\ln \stackrel{\mathrm{‾}}{L}-\mathit{\mu }{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}},\stackrel{\mathrm{‾}}{L}>\mathrm{0},\end{array}$$

Simulations and analyses were conducted corresponding to the $\mathrm{8}\times \mathrm{6}\times \mathrm{6}\times \mathrm{18}\times \mathrm{13}$ (684 400) variable combinations for the DFNs, with each case independently simulated 20 times. The number of variations for each variable are listed in Table 5.

Table 5 Parameters of DFNs in Sect. 3.3.

3.4 Derivation of percolation threshold equation for DFN with lognormal fracture lengths

The multivariable fitting process for the DFN with lognormal fracture lengths was analogous to that described in Sect. 3.2. Initially, the hypothetical relationship between n_t and $\stackrel{\mathrm{‾}}{L}$ was fitted (Eq. 18). Subsequently, by analysing the relationship between the parameters in the hypothetical model, new variables ν, Δ, and κ were sequentially incorporated. The expressions for the newly added variables were then used to replace the parameters in the hypothetical model. Ultimately, this yields the fitted relationship between n_t and $\left(\stackrel{\mathrm{‾}}{L},\mathit{\nu },\mathrm{\Delta },\mathit{\kappa }\right)$, with the fitting process detailed in Eqs. (18)–(21), where x=cos Δ.

$$\begin{array}{}\text{(18)}& {\displaystyle }{n}_{\mathrm{t}}=f\left(\stackrel{\mathrm{‾}}{L}\right)=a{\stackrel{\mathrm{‾}}{L}}^b+c,\text{(19)}& {\displaystyle }{n}_{\mathrm{t}}=f\left(\stackrel{\mathrm{‾}}{L},\mathit{\nu }\right)=a_{\mathrm{1}}{\mathrm{e}}^{a_{\mathrm{2}}\mathit{\nu }}{\stackrel{\mathrm{‾}}{L}}^{b_{\mathrm{1}}}+c_{\mathrm{1}}{\mathit{\nu }}^{c_{\mathrm{2}}},\text{(20)}& {\displaystyle }\begin{array}{rl}{n}_{\mathrm{t}}& =f\left(\stackrel{\mathrm{‾}}{L},\mathit{\nu },\mathrm{\Delta }\right)=\left(a_{\mathrm{11}}x+a_{\mathrm{12}}\right){\mathrm{e}}^{a_{\mathrm{22}}\mathit{\nu }}{\stackrel{\mathrm{‾}}{L}}^{b_{\mathrm{11}}x+b_{\mathrm{12}}}\\ & +c_{\mathrm{11}}{\mathit{\nu }}^{c_{\mathrm{12}}x+c_{\mathrm{13}}},\end{array}\text{(21)}& {\displaystyle }\begin{array}{rl}{n}_{\mathrm{t}}& =f\left(\stackrel{\mathrm{‾}}{L},\mathit{\nu },\mathrm{\Delta },\mathit{\kappa }\right)=\left[(d_{\mathrm{1}}{\mathit{\kappa }}^{\mathrm{2}}+d_{\mathrm{2}}\mathit{\kappa }+d_{\mathrm{3}}\right)x\\ & +(d_{\mathrm{4}}{\mathit{\kappa }}^{\mathrm{2}}+d_{\mathrm{5}}\mathit{\kappa }+d_{\mathrm{6}})]{\mathrm{e}}^{\left(d_{\mathrm{7}}\frac{\mathit{\kappa }}{\stackrel{\mathrm{‾}}{L}}+d_{\mathrm{8}}\right)\mathit{\kappa }\mathit{\nu }}{\stackrel{\mathrm{‾}}{L}}^{d_{\mathrm{9}}{x}_{+}^{\mathrm{2}}{d}_{\mathrm{10}}\mathit{\kappa }+d_{\mathrm{11}}}\\ & +d_{\mathrm{12}}{\mathit{\nu }}^{d_{\mathrm{13}}{x}^{\mathrm{2}}+d_{\mathrm{14}}},\end{array}\end{array}$$

Similarly, the parameters are used as inputs for the final multivariable least squares optimisation based on the Levenberg-Marquardt algorithm, utilising all simulation results. The final optimised values of the parameters in Eq. (21) are presented in Table 6.

Table 6 Parameters of the percolation equation in terms of fracture properties.

To visualise the relationships in Eq. (21), several surfaces of n_t vs. $(\stackrel{\mathrm{‾}}{L},\mathit{\nu },\mathrm{\Delta },\mathit{\kappa })$ corresponding to different values of ν (2, 4, and 12) and κ (2, 7, and 10) are presented in Fig. 13. Generally, higher Δ and lower $\stackrel{\mathrm{‾}}{L}$ result in increased percolation threshold n_t.

Figure 13 Extracted surfaces from the final percolation threshold equation. (a) $\mathit{\nu }=\mathrm{4},\mathit{\kappa }=\mathrm{2}$; (b) $\mathit{\nu }=\mathrm{4},\mathit{\kappa }=\mathrm{7}$; (c) $\mathit{\nu }=\mathrm{2},\mathit{\kappa }=\mathrm{10}$; (d) $\mathit{\nu }=\mathrm{12},\mathit{\kappa }=\mathrm{10}$.

Equation (21) is derived for a dimensionless unit square, with its results expected to be applicable to regions of varying scales (y×y). For these scenarios, the percolation threshold ${\stackrel{\mathrm{‾}}{n}}_t$ will be adjusted using a scaling modification from Eq. (21), resulting in Eq. (22). If the number of fractures N is more than ${\stackrel{\mathrm{‾}}{n}}_{\mathrm{ts}}$, the probability of achieving percolation will be high.

$$\begin{array}{}\text{(22)}& n_{ts}=n_{\mathrm{t}}\times y_n,\end{array}$$

where y_n is the scaling correction factor, $y_n={\log }_{\mathrm{10}}(\frac{\mathrm{1}}{y}+\mathrm{20})$.

4 Validation of the derived percolation threshold equation

4.1 Percolation analysis of fracture networks not used for deriving the threshold equation

1. DFN with exponential fracture length. To test the performance of the derived zero percolation thresholds (Eq. 14), additional DFN models with different parameters ($n=\mathrm{100},\mathrm{200}$, Δ=3°, 63°, $\mathit{\kappa }=\mathrm{5},\mathrm{40},\mathrm{50})$ at different scales (2 m×2 m, 60 m×60 m, 300 m×300 m, 900 m×900 m, 1100 m×1100 m, 1200 m×1200 m, 1300 m×1300 m) were generated for percolation analysis. The percolation thresholds obtained from Eq. (14) and numerical simulation are shown in Fig. 14. The close agreement between the predicted thresholds and the simulation results demonstrates that the derived relationships (Eq. 14) performed extremely well for predicting of zero percolation thresholds of DFNs with different parameters at different study scales.

2. DFN with lognormal fracture length. To test the derived zero percolation thresholds (Eq. 22), additional DFN models with different parameters ($\stackrel{\mathrm{‾}}{L}=\mathrm{0.150.25}$, $\mathit{\nu }=\mathrm{5},\mathrm{7}$, $\mathrm{\Delta }=\mathrm{26}\mathit{°},\mathrm{56}\mathit{°}$, $\mathit{\kappa }=\mathrm{5},\mathrm{9})$ at different scales (2 m×2 m, 60 m×60 m, 300 m×300 m) were generated for percolation analysis. The percolation thresholds obtained from Eq. (22) and numerical simulations were depicted in Fig. 15. The remarkable concurrence between the predicted thresholds and the simulation outcomes underscores the robust performance of the derived relationships (Eq. 22) in accurately forecasting percolation thresholds across diverse parameter configurations and study scales for DFNs.

Figure 14 Validation of Eq. (14) in study areas on different scales.

Figure 15 Validation of Eq. (22) in study areas on different scales.

4.2 Comparison with analytical solutions of simple fracture networks

The equations derived in Sect. 3.2 are for stochastic fracture networks with fracture lengths that follow an exponential distribution and fracture orientations that follow a von-Mises distribution. Because of the complexity of such fracture networks, there is no analytical solution for the corresponding percolation threshold. However, for simple fracture networks, where both fracture location and orientation follow completely random distributions and the fracture length is identical, the analytical solution for the percolation threshold of fracture length is (Balberg et al., 1984; Berkowitz, 1995):

$$\begin{array}{}\text{(23)}& L_C=\mathrm{4.2}/\sqrt{\mathit{\pi }\mathit{\rho }},\end{array}$$

where ρ is the fracture density (= P20) calculated as $\mathit{\rho }=n/y^{\mathrm{2}}$ and y² is the area of the study region. For the case of varying fracture length, the corresponding threshold is:

$$\begin{array}{}\text{(24)}& {\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=\sqrt{L_C^{\mathrm{2}}-{\mathit{\sigma }}^{\mathrm{2}}},\end{array}$$

where σ² is the variance of fracture length distribution. If the length follows an exponential distribution with parameter λ, Eq. (24) becomes (Berkowitz, 1995). Therefore, this section utilises fracture networks characterized by an exponential distribution of fracture lengths as a case study to compare the derived thresholds with analytical solutions.

$$\begin{array}{}\text{(25)}& {\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=L_C/\sqrt{\mathrm{2}}=\mathrm{4.2}/\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\pi }n}{y^{\mathrm{2}}}}}\approx \mathrm{1.676}y{n}^{-\mathrm{0.5}},\end{array}$$

This is a special case covered by the relationships derived in this work by setting κ=0 for a completely random distribution of fracture orientation and ignoring cos Δ as it is now irrelevant. Eq. (14) then becomes:

$$\begin{array}{}\text{(26)}& {\stackrel{\mathrm{‾}}{L}}_{\mathrm{ts}}=y{a}_{\mathrm{32}}{n}^{b_{\mathrm{33}}},\end{array}$$

where a₃₂=0.9823, $b_{\mathrm{33}}=-\mathrm{0.4730}$ The equation can be further simplified to:

$$\begin{array}{}\text{(27)}& {\stackrel{\mathrm{‾}}{L}}_{\mathrm{ts}}=\mathrm{0.9823}y{n}^{-\mathrm{0.4730}},\end{array}$$

which should be compared to Eq. (25). Note that the difference between the two equations is due to the different probabilities used to derive the percolation threshold. The theoretical solution is for a percolation probability of 50 % while the derived relationship is for a percolation probability of 0 %, as discussed above, and therefore it should be smaller.

There is a striking similarity between the analytical solution for p50 and the solution we derived for p0. The reasons why these two equations are so similar and what the factor of two in the first coefficient represents are important. These will be discussed in future work instead of here since it is not the focus of this work.

To compare the solutions, fracture networks with $n=\mathrm{100},\mathrm{200},\mathrm{300},\mathrm{400}$ and 500 in a square of 75 m×75 m are used for the simulations. Percolation thresholds corresponding to p, p50 and p100 are calculated by Monte Carlo simulations and the results are shown in Fig. 16. As demonstrated, the analytical solution is close to the p50 percolation threshold with an average absolute difference of 4.9 %. On the other hand, the solution based on the derived equation is close to that of p with an average absolute difference of 11.4 %.

Figure 16 Comparison of percolation thresholds determined by the analytical solution, the derived equation, and numerical simulations. The points are the averages of 20 groups of MC simulations and the error bars are three times the corresponding standard deviations.

4.3 Percolation analysis of real fracture networks using the derived equation

Two real fracture networks, as shown in Figs. 17a and 18a, are used to demonstrate further the application of the derived percolation threshold equations. Figure 17a shows a set of fractures traces on a rock outcrop taken from Wilson (Wilson, 2001). Figure 18a shows fracture traces in the deformation bands on the Valley of Fire State Park, Nevada (Barton and Hsieh, 1989). Mid-points of the fractures are used to represent the fracture locations, as shown in Figs. 17b and 18b, respectively. They are all considered to follow approximately the Poisson distribution. The numbers of fractures, n, in the two systems are 35 and 186. Clearly there is one dominant direction of fracture orientations in these systems, as illustrated in the rose diagrams shown in Figs. 17c and 18c. The orientation dispersion parameters (κ) were calculated to be 145.18 and 25.68. For fracture length, the histogram in Fig. 17d indicates an approximately exponential distribution for the first fracture set. For the second and third sets, the histograms (Fig. 18d) suggest lognormal distributions. The average fracture lengths $\stackrel{\mathrm{‾}}{L}$ are 0.249 m, 0.282 mm, 0.9544 m and the side lengths y of study areas are 2.5 m, 8 mm, 7 m, respectively.

Figure 17 A set of fractures from an outcrop (Wilson, 2001) and the network properties. (a) Fracture network from an outcrop; (b) Fracture centres; (c) Fracture orientation; (d) Fracture length.

Using Eqs. (21) and (22) with the parameters $(\stackrel{\mathrm{‾}}{L},\mathit{\nu },\mathrm{\Delta },\mathit{\kappa })$ listed in Table 7, the calculated percolation thresholds (y_n×n_t) in the horizontal and vertical directions can be calculated. For Fig. 17a, the calculated percolation thresholds are 126 (horizontal) and 94 (vertical). The threshold in the horizontal direction is much greater than that in the vertical direction due to the fact that fractures are mainly vertical in this case. The fracture numbers are all less than these two thresholds, hence the fracture network is not percolated in both directions. This conclusion can easily be confirmed in this case by visual inspection of the fracture system displayed in Fig. 17a.

Table 7 Parameters of the real fracture network in Fig. 17a and percolation assessment.

For the fracture set (Fig. 18a), the average fracture length is 1.19 m. The horizontal percolation threshold $y\times {\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}$ is 1.34 m and the vertical threshold is 0.69 m. The average fracture length in this case is greater than the vertical threshold and therefore the fracture network is percolated vertically but not horizontally. On close inspection of Fig. 18a, there is a cluster of fractures (marked in red) connecting the top and bottom sides of the study region.

Table 8 Parameters of the real fracture network in Fig. 18a and percolation assessment.

Figure 18 Fracture traces of deformation bands in the Valley of Fire State Park, Nevada (Barton and Hsieh, 1989) and the corresponding network properties. (a) Fracture network from an outcrop; (b) Fracture centres; (c) Fracture orientation; (d) Fracture length.

5 Discussions

While the simplified relationships proposed in this study provide an efficient means of estimating the zero-percolation threshold directly from fracture network parameters, several limitations should be acknowledged. First, the present formulation is restricted to a two-dimensional (2D) domain with a single fracture set of similar orientation. This abstraction neglects the inherently three-dimensional (3D) nature of natural fracture systems, which typically comprise multiple intersecting sets formed under polyphase tectonic regimes. Extending the methodology to accommodate 3D configurations with multiple fracture sets will be a critical next step, although the associated analytical derivations are considerably more complex.

Second, the current approach relies on specific statistical assumptions for fracture attributes (e.g., Poisson-distributed locations, exponential or lognormal length distributions, and von Mises orientation distributions). While these choices are consistent with many prior DFN studies, natural fracture systems often exhibit spatial correlations, clustering, or heavy-tailed distributions that may deviate from these idealizations. Future work will explore non-parametric estimation methods to mitigate dependency on pre-defined distribution forms and to better capture geological variability.

Third, the present formulation addresses only the geometric aspects of network connectivity. In real-world applications, hydraulic connectivity depends on both geometry and fracture transmissivity, which is influenced by aperture, roughness, and in-situ stress. Incorporating coupled flow–mechanical models would enable the framework to more directly predict hydraulic properties such as permeability.

A particularly promising avenue lies in integrating the derived functional forms of percolation thresholds into physics-informed machine learning workflows. The explicit parameter, threshold relationships, established here can serve as strong physical constraints within emerging architectures, such as Kolmogorov–Arnold Networks (KAN), to enhance generalizability and reduce training data requirements. Embedding these analytical insights into AI frameworks is expected to significantly improve predictive efficiency for complex, data-limited geological systems.

While the proposed 2D, single-set formulation constitutes a foundational step, its extension to multi-set, 3D, and distribution-flexible fracture systems, together with integration into advanced machine learning paradigms, offers a clear and impactful trajectory for future research.

6 Conclusions

Percolation analysis of fracture networks is important for many applications, including oil and gas recovery, geothermal energy exploitation, hydrology, and groundwater protection in radioactive waste storage. In this paper, we focus on the percolation threshold relevant to rock impermeability, which is critically important for the safe underground storage of waste and energy materials.

Our approach to the calculation of the percolation threshold makes direct use of the characteristic parameters of 2D fracture network, in particular the number of fractures n, the fracture size (length) $\stackrel{\mathrm{‾}}{L}$ and the fracture orientation Δ. This differs from the simplified approaches of using indirect characteristic parameters (e.g., fractal dimension), which could produce misleading results because fracture orientation is not considered. The assessment of fracture networks in this research was made under the following assumptions: (1). the centre points of fractures are randomly and independently distributed in space; (2). the lengths of fractures follow an exponential distribution; and (3). the orientations of fractures follow a von-Mises distribution, the parameters of which are the mean orientation μ and the concentration parameter κ. The relationship between fracture network parameters and the corresponding percolation threshold is obtained from a large number of simulations. A non-linear multivariate fitting process was used to derive the final prediction equation for the percolation threshold in the form of ${\stackrel{\mathrm{‾}}{L}}_{\mathrm{t}}=f(n,\cos \mathrm{\Delta },\mathit{\kappa })$. The derived equation provides a reliable relationship and an efficient way to estimate the connectivity and percolation state of a fracture network based directly on its parameters. The relationship was cross-validated using a published analytical solution and was further applied to three real fracture networks. The results demonstrate that the derived relationship can be used for fracture networks at different scales using a rescaling coefficient and can also be used for the assessment of percolation in different directions. The derived relationship is a useful extension for rock impermeability evaluation (zero probability percolation p0), compared with the commonly used percolation assessment based on excluded volume, which corresponds only to the occurrence of percolation on average (i.e., 50 % probability percolation, p50). Additionally, this work also studies fracture network models with log-normally distributed fracture lengths and derives zero percolation formulas, reaching conclusions similar to those mentioned above.

Owing to the complexity of multiple fracture sets, the present study is confined to a 2D fracture network with a single set. Future work will extend the framework to both 2D and 3D systems with multiple fracture sets, incorporating machine learning techniques to address the increased complexity.

Appendix A: Description of parameters

The parameter descriptions are shown in Table A1.

Table A1 Description of parameters.