Analysing the Connectivity of Fracture Networks Using Natural Fracture Characteristics in the Khairi Murat Range, Potwar Region, Northern Pakistan

1. Introduction

Fractures are ubiquitous in rocks and have long been recognized as an essential property of rock mass due to their numerous applications in various fields. (Alghalandis, 2017). For example, fractures, being as essential discontinuities significantly influence the rock quality in engineering geological applications (Andersson et al., 2002); the strength of rock masses (Hoek & Martin, 2014; Cao et al., 2018; Lin et al., 2019); slope stability (Hadjigeorgiou & Grenon, 2005); and civil structures on the earth’s surface such as dams and foundations. (Staub & Outters, 2002; Alfven, 2015; Lin et al., 2019). The fracture network and its potential connectivity factors are indispensable for the safe building and development of underground repositories (Faoro et al., 2016), as well as the storage (Ye et al., 2021) or disposal of hazardous wastes. (Follin et al., 2006; Karyounis & Jenny, 2011; Reeves et al., 2012; Alghalandis et al., 2013; Chang et al., 2017). Rock fracture networks are important variables that influence fluid flow through rock masses in the subsurface. (Dverstorp, 1991; Pan et al., 2010). For groundwater transportation and reservoir hydrocarbon potential, rock fractures and interconnectivity of fracture networks are extremely important (Dejam et al., 2018; Li et al., 2017; Berkowitz, 2002; Koyama et al., 2009; Acharya et al., 2014). This research employs a modified version of the methods outlined by Ghosh & Mitra (2009) to characterize the rock fracture properties in terms of connectivity of fracture networks within rock masses as a medium for fluid transport.

Rock fractures are analyzed by the characterization of multiple fracture parameters such as orientation, spacing, aperture, density or frequency, length, roughness and so on. Several researchers extensively use these parameters for their particular objectives, and are well documented in the literature. Fracture orientation, density, and length play an important role in determining the connectivity of fracture networks and predicting flow pathways because the connectivity is a measure of the interconnectedness of fractures and fracture networks in a system domain (Pedretti & Bianchi, 2019; Ye et al., 2021). Connectivity of the fracture networks in fractured rock generally depends on three parameters: the degree of connectedness of individual fractures (based on the relative dispersion of the orientation of fracture sets); the number of fractures per unit area called density in a 2D framework, and fracture trace length (Wu & Pollard, 2002; Leckenby et al., 2005; Sanderson & Nixon, 2018). Several research studies (Hadgu et al., 2017; Somr et al., 2017) have been carried out in an attempt to gain a better understanding of how fracture attributes and their statistical properties contribute to the connectedness of fracture networks (Zhang & Sanderson, 2002; Lee et al., 2009; Zhang et al., 2013; Hadgu et al., 2017; Jiao et al., 2018).

Numerous published studies, including the one referenced herein, have proposed various conceptual models such as the discrete Fracture Network (DFN) and the Equivalent continuum Model (ECM) to address the complexities of fluid flow and the diversity of fracture connectivity patterns. These models incorporate different methodologies for simulating fractures, estimating the probabilistic distribution of fracture properties, and generating the resulting fracture networks. Fracture attributes datasets are commonly assumed, synthesized, or integrated using computational techniques to meet model input requirements. The resulting data are then extrapolated across the targeted natural fracture host domain. A wide range of conceptual and numerical models, along with specialized computer codes and software, are being developed rapidly. In recent years, several new modeling and simulation methodologies have been introduced into the literature to analyze systems and their associated geometric parameters (Xu et al., 2006; Sausse et al., 2008; Molebatsi et al., 2009; Tavani et al., 2010; Xu & Dowd, 2010; Ghaffari et al., 2011; Alghalandis et.al., 2013; Le Goc et.al., 2014; Hyman et.al., 2015; Maillot et.al., 2016; Alghalandis, 2017; Lavoine et.al., 2020; Nyberg et al., 2018). Although exponential advancements in fracture modelling and simulation offer significant advantages and pave the way for innovative research and modern methodologies, the importance of field investigations at the outcrop scale remains indispensable. For novice researchers, the process of collecting geological data directly from the field is often more appealing than the complex and sometimes overwhelming task of selecting optimal models or developing new software and computational tools for simulation and analysis. Field work not only provides tangible insights into geological features but also fosters a deeper, more intuitive understanding that purely computational approaches may lack.

Therefore, to address this problem, we used an analytical and investigative approach based on field data supplemented with numerical values by mathematical equations and statistical methods. Fracture network characterization and fracture connectivity analysis are based on the interpretation of integrated fracture data obtained from field measurements, geo-structural maps and fracture trace maps developed from field photographs. The objective of this paper is to evaluate and identify the role of fracture properties in the characterization of fracture connectivity using numerical and statistical distribution analyses. The area of Khairi Murat (Fig. 1) was selected as a case study for this purpose because it contains excellent outcrop exposures with remarkable imprints of structural features, such as fractures, faults, and folds.

2. Geology and Structure

2.1. Geology of the Area

The Khairi Murat Range (KMR) is located about 40 kilometers southwest of Islamabad, the capital city of Pakistan, between the northern limb of the Soan syncline in the south and the Main Boundary Thrust (MBT) in the north. The area is part of the northeastern edge of the Potwar Sub Basin and contains a stratigraphic succession from the Early Eocene to the Middle Pliocene (Kazmi & Jan, 1997; Kazmi & Abbasi, 2008) (Fig. 2). The Early Eocene assemblage of Charrat Group comprises Margallah Hill Limestone, Chorgali Formation and Kuldana Formation which collectively preserve a record of shallow marine, marine to continental paleo-environment (Khan et al., 1991; Mujtaba, 2001; Benchilla et al., 2002; Jurgan & Abbas, 2006; Khan et al., 2006; Zawar et al., 2017; Awais et al., 2020). Due to a period of non-deposition in the Potwar Fold Belt (PFB), the Early Eocene to Oligocene stratigraphic sequence is missing (Tab. 1). The Miocene Rawalpindi Group includes Murree Formation and Kamlial formations, is extensively exposed in the study area. The lithological units of the Rawalpindi Group consist of thin to thick bedded sandstone, cross-bedded siltstone, and mudstone, with remarkable fluvial depositional imprints. The Chinji, Nagri, and Dhok Pathan Formations of the Siwalik Group designate the Mio-Pliocene molasses sequence of sediments that are exposed in the study area (Shah, 1977; Ali et al., 2023). Outcrops of these stratigraphic sequences preserve a significant record of tectonic events in the form of geological structures.

2.2. Structure of the Area

Structurally, the Potwar Plateau is divided into two zones: the North Potwar Deformed Zone (NPDZ) and the South Potwar Platform Zone (SPPZ) (Fig. 3a). The study area is located in the North Potwar Deformed Zone (NPDZ). This zone is a product of the Himalayan orogenic deformation system and experiences complex deformation including imbricate thrusts, ramp-flat duplexes, and triangle zones developed during the Neogene (Jaswal et al., 1997; Jadoon & Frisch, 1997; Jadoon et al., 1997, 2008, 2015). In this zone, the sedimentary sequences are imbricated by southward thrusting and fault-related folding, with a minimum crustal shortening rate of 22mm per year (Moghal et al., 2007; Jaswal et al., 1997; Iqbal & Bannert, 1998). The southern boundary of NPDZ is marked by the Khairi Murat Range, which has been tectonically uplifted along the south-verging Khairi Murat Thrust. This high-angle fault, striking NE-SW, runs approximately perpendicular to the tectonic transport direction. The fault is clearly evident where competent Eocene strata overlie the Kuldana and Murree Formations (Ali et al., 2023).In the footwall of the Khairi Murat Thrust to the south, the Siwaliks have been thrust over by the Dhurnal fault, a north-facing roof thrust. A triangle zone exists between the Khairi Murat Thrust and the passive back thrust (Dhurnal fault) (Jadoon et al., 1997; Jaswal et al., 1997; Jaswal, 1999; Kamran & Siddiqi, 2011; Jadoon & Frisch, 1997; Jadoon et al., 2015). The Khairi Murat Anticline is a significant fold in the study area. Its core is composed of Margalla Hill limestone, while its northern limb comprises the Chorgali, Kuldana and Murree Formations. The southern limb is faulted (Jadoon & Frisch, 1997) (Fig. 3). The study area also contains several small-scale strike-slip offsets, such as the Gali Jagir-Chauntra fault (Fig. 3a and 3b) and records approximately 55km of crustal shortening, indicating the complexity of tectonic deformation in this region (Baker et al., 1988; Pennock et al., 1989).

3. Data Acquisition Techniques and Analysis Methods

3.1. Data Collection from Outcrops

A rock fracture sampling operation was conducted at 107 sites within a 35 km × 2.5 km area of the Khairi Murat Range. All possible lithological and structural domains were surveyed. A hand-held GPS device is used to record the positions of the sampling sites so that the fracture maps can be georeferenced. Approximately 3,224 individual fractures were mapped across 107 sites in 17 data collection centers (Fig. 3b shows the locations of these centers). A series of maps was produced by taking photographs directly from outcrops, focusing on the fractures in the field in sedimentary strata. Two techniques were employed to obtain representative sampling in the study area.

3.1.1. Circular Scanline Technique

The circular scanline method (Baecher, 1980; Mauldon, 1998; Mauldon et al., 2001; Rohrbaugh et al., 2002), also referred to as the circle inventory window (Davis & Reynolds, 2011), is a significant and effective tool for fracture data collection. It offers a wide range of applications and provides time-efficient estimates of trace density and mean trace length (Mauldon et al., 1999). A circular scanline is simply a circle drawn on the surface of a rock (Fig. 4b). Based on insights and guidance from Zeeb et al. (2001), Mauldon et al. (2002), Rohrbaugh et al. (2002), IBN (2000), and Kamali et al. (2016), and following the arguments of Davis and Reynolds (2011), a fixed-radius circle (sample size) is preferred for measuring fracture properties.

A fracture trace map is created by physically tracing visual fracture signatures directly on outcrop faces. This method minimizes sampling bias by reducing complications associated with the random distribution of fracture lengths and anisotropy in orientations, thereby enhancing practical applicability. The estimator equation developed by Davis and Reynolds (2011) is used to determine fracture density from outcrops and photographs (i.e., photographs with known focal length, height, and scale) (Fig. 4c & d).

3.2. Data Acquisition in the Lab

3.2.1. Photograph Trace Maps

AutoCAD software provides an effective working environment for the characterization and quantification of natural fracture features. All data were digitized and analyzed using this software. Based on shared geometrical and physical properties, such as orientations (with fixed dispersion intervals), dip magnitude, and dip direction, individual fractures were grouped and presented as rose diagrams (Fig. 5)

4. Analysis Methods

4.1. Method adopted from Davis & Reynolds (1997)

The density of fractures at a sampling station is defined as the number of fractures per unit area, as measured using the circle inventory method (Dershowitz and Herda, 1992; Ghosh and Daemen, 1993). Fracture density depends on three primary rock mass properties: the mechanical property of rock unit, bed thickness and proximity to structural domain (i.e., intensity of applied stress) (Goodman & Shi, 1985; Florez-Nino et.al., 2005). In this study, fracture density is quantified as the total number of fractures (regardless of fracture sets) occurring within a unit-radius circle, in accordance with the circle inventory method. The value is mathematically derived using the equation provided below (Davis & Reynolds, 2011):

$$Pf={\displaystyle \raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{${\pi r}^2$}\right.}$$

(1)

Where,

$Pf$= Fracture density

$L$= Cumulative length of all fractures

$r$= radius of the inventory circle

The measure of density is expressed as length per unit area (e.g., ft/ft², cm/cm², m/m², km/km²) or in reciprocal form. (i.e., ft^-1, cm^-1, m^-1 and km^-1). Density values (Tab. 3) play an important role in establishing the connectivity of fracture networks.

4.2. Method adopted from Ghosh & Mitra (2009)

The interconnectedness of fractures is determined by fracture characteristics, as described in Section 1. Fracture length, density, and orientation are three major fracture characteristics that influence the connectivity of fracture networks. Interconnected fractures form clusters, which are then linked to nearby clusters by one or more longer fractures. These linking fractures do not allow the fracture networks to remain isolated and are termed as “backbone” fractures (Alghalandis et al., 2013; Alghalandis, 2014; Zhu et al., 2021), previously termed “skeleton” fractures by Priest (1993) and Liang (2016). The interconnection of individual fractures into clusters, along with the linkage of these clusters via backbone fractures, results in fully connected fracture networks. Such networks represent the overall connectivity of the rock fracture system (Stauffer, 1985; Berkowitz & Balberg, 1993; Berkowitz, 1995; Ghosh & Mitra, 2009).

To estimate the degree of fracture connectedness within the network, a cluster analysis technique adopted from Ghosh and Mitra (2009) is used.

The mean connectivity and general characteristics of fracture networks in an area can be characterized by the total length of connected fractures (FCL) or the area of connected fracture networks (FCA). Fracture cluster length (FCL) can be expressed as:

$$\mathrm{FCL}₌{\displaystyle \frac{Totallengthofconnectedfractures}{Totalsamplearea}}$$

(2)

Whereas fractional connected area (FCA) is:

$$\mathrm{FCA}₌{\displaystyle \frac{Totalsurfaceareaofrockconnectedbyfractures}{Totalsamplearea}}$$

(3)

The FCA is generally considered the most suitable way of evaluating the effectiveness of fracture networks in facilitating fluid flow (Odling, 1992, 1995, 1997; Ghosh & Mitra, 2009). To estimate the FCA, high-resolution photographs of outcrops exhibiting the best fracture exposure were taken at 67 sites. The analysis was performed using the open-source software Open Plot (Tavani et al., 2011).

The connectivity of the fracture network is often evaluated using percolation theory, which is based on the interconnectivity of fractures and fracture network density, fracture length, and orientation. (Molebatsi et al., 2009). Complex regional and local stress perturbations generate multiple fractures and fracture sets, including extension and hybrid shear fractures with wide orientation dispersion (Hancock et al., 1987). Based on the statistical Fisher distribution (Mardla, 1972), eight fracture sets were identified. Orientation measurements taken from field outcrops were grouped into fracture sets using the techniques outlined by Michelena et al. (2013) and Zangerl et al. (2022), in conjunction with structural geological criteria derived from field observations (Ozkaya et al., 2011; Jiang et al., 2016). Fracture orientations (defined by the strike and dip trends) are used to describe fracture clustering behavior and its effect on network connectivity.

5. Analysis and Results

The region of interest includes complex structural features such as thrusts, faults, folds, deformation bands, and strike-slip components associated with episodic deformation. A significant number of fracture populations related to these deformational features are present in the study area (Dasti et al., 2018). Based on the explanation of Grodner et al. (2021) and assuming a constant fracture aperture, it is observed that nearly all fractures measured at sampling stations are open and fluid-conducive fractures. The analytical methods used in this study, when combined with statistical analysis, become twofold helpful tools for investigating fracture network connectivity.

5.1. Fracture Orientation Analysis

Eight fracture sets are identified at 67 surveyed locations in the study area using circular statistics and the methodology described by Ozkaya (2011). Based on their physical characteristics, the entire fractures population was classified into sets considering a) similarities among strike and dip, and b) the mode of deformation with angular dispersion of 20° for strike and 10° for dip of each set as shown in the rose diagram (Fig. 5). Table 2 shows the mean strike and dip of these eight fracture sets: N4°W/87°(vertical), E-W/49°NW, N55°E/50°SE, N30°E/90(vertical), N45°E/49°NW, N56°W/45°NE, N43°W/90(vertical), and N30°W/88°(vertical). These fractures developed in response to their proximity to structural elements and reflect the episodic occurrence of various tectonic events. A broad range of angular scatter in fracture orientation increases the likelihood of individual fractures intersecting to form connected networks. Out of the eight fracture sets, nearly two to three are observed at each measurement station.

Fracture Set-1, which runs from N4°W to N20°E, is the dominant set in the study area, followed by Set-8, the second-largest set with a narrow range of strike angles between N25°W to N35°W. Several fractures in Set-1 strike between N20°W and N20°E, including a north-south (N-S) orientation, with slight distortion in dip and strike angles. Similarly, in Set-2, a number of individual fractures occur along E-W or at orientations of N85°W and N85°E. Set-3 strikes in NE-SW, and Set-4 in ENE-WSW. Set-5 follows the NNW-SSE directions, while Set-6 follows NW-SE. Set-7 trends N40°-45°W. In the study area, there are two groups of orthogonal fracture systems (Set-1 & Set-2 and Set-4 & Set-7) and one group of conjugate fracture systems (Set-3 & Set-5). The fractures in Set-3 and Set-5 are in the same quadrant but dip in opposite directions, i.e., 50°SE and 49°NW, respectively (Tab. 2). The rose diagram depicts the manual grouping of the total measured population of fractures (Fig. 5). The Box-and-Whisker plot of the orientation data (Fig. 6a) shows that the fracture population of Set-1 is more scattered than those of the other sets of fractures. The frequency distribution of each group is presented in Fig. 6b. The statistical characterization of fracture density and clustering in terms of fractional connected area (FCA) provides a basis for assessing fracture network connectivity.

5.2. Fracture Density Analysis

To avoid sampling bias, as described in Rohrbaugh et al. (2002), a uniform circle (with equal radius, i.e., 1m) was drawn on the bedding surface of exposed outcrops throughout the sampling stations. The entire region was divided into five zones, as shown in the fracture density map (Fig. 7). In Zone I, the fracture density varies from 0.5 to 5 fractures per meter, with a mean density of 4.63. Zone II has a density of 6-10 fractures per meter with a mean density of 8.57. Similarly, Zones III and IV have concentrations ranging from 11 to 16 fractures per meter and 16 to 20 fractures per meter, with mean densities of 13.37 and 18.61, respectively. Zone V has a higher density and varies from 20 to 27 fractures per meter, with a mean density of 22.81, as summarized in Table 3. A visual inspection of the density map (Fig. 7) and the superimposition of density values on the structural map of the study area demonstrate that density has a significant correlation with structural elements and deformation intensity. Table 4 shows the statistical properties of cumulative fracture density across the entire region.

5.3. Fracture Length Analysis

The fracture length distribution for each field-scale fracture group has been examined and quantified using values obtained from both the circle inventory method and fracture trace maps derived from field photographs. Cumulatively, the eight fracture sets have a mean fracture length of 11m with a minimum of 3m and a maximum of 27m. As illustrated in Fig. 9 (frequency distribution), the majority of fractures remain in the length bracket of 6m to 14m. Shorter fractures are more abundant, while longer fractures are scarce and occur only sparsely in the study area. Similar to fracture density, fracture length is also strongly influenced by structural position and lithological contrasts.

5.4. Connectivity and Cluster Analysis

The 2-D fracture trace maps were generated by manually tracing fractures on photographs. The connected fractures were characterized and grouped into clusters and connected areas. Fracture trace maps prepared from a satellite imagery were also used to analyze scale-dependent variation in the estimation of fractional connected area (FCA). Statistics for FCA, referred to as “connectivity” in this paper, are presented in Table 5. As shown in Table 5, the connectivity of fractures in the study area ranges from a minimum of 0.070 to a maximum of 0.820, with a mean value of 0.404. These estimated connectivity values are considered relative, as absolute values cannot be directly measured from the outcrops.

The number of fracture sets greatly influences the FCA and determines the size of clusters in a given area. The distribution of fracture orientations facilitates the formation of cluster size and affects connectivity. Fracture density is also closely associated with the formation of cluster size and the number of fractures in the system. Furthermore, while the scale and area of observation have minimal effect on properties of individual fractures (e.g., orientation, density, length), they greatly influence properties related to the fracture systems, such as network geometry and connectivity. Figure 10 and Figure 11 represent two sets of fractures. An analysis of these figures clearly shows that clusters in Fig. 10 are bigger than those in Fig. 11. In Fig. 10, the proportional connected area (FCA) of fracture sets is calculated to be 0.30, whereas in Fig. 11, it is 0.125, which is smaller than in Fig. 10. Similarly, the fracture density in Fig. 10 is 10 fractures per meter, which is greater than the density in Fig. 11, which is 4 fractures per meter. Figure 12 shows three fracture sets on the same scale as Figure 10 and Figure 11, but the observational region varies. The interconnected area (FCA) goes to 0.71 in this figure (Fig. 12), and the estimated density rises to 18 fractures per meter. According to Fig. 13, the fracture-connected area is 0.74, which is similar to the area shown on a smaller scale (Fig. 12), and the density is dependent on the number of fractures (7 fractures per meter) and the dispersion of fracture orientation. Longer fractures dominate as the observational region and scale increase. In Fig. 14, it is clear that the longer fractures play a key role in network connectivity (cluster interconnection) and serve as “backbones” in the system’s interconnection process to some extent, but are limited to indirect control over FCA, density, and cluster size. The role of longer fractures in connectivity and density variations shown in these examples is scale-dependent. The data collected directly from the outcrops are presented in Fig. 15 and demonstrate a non-linear relationship between the density of fractures and the length of individual fractures.

As shown in Fig. 16, fracture connectivity is plotted against fracture density (per meter). It is evident that connectivity is directly proportional to density and can be adequately characterized by a linear function with an R² close to 0.97. Figure 17 demonstrates that the length data is dominated by short-length fractures and that there is only a weak correlation between connectivity and fracture length. The connectivity vs length plot (Fig. 17) shows that connectivity decreases as the length of individual fractures increases. In the length range of 0-100cm, connectivity declines, and as the fracture length exceeds 100cm, the impact on connectivity diminishes.

6. Discussion

As Renshaw (2000) explains, the statistical properties of fracture parameters, on which fracture networks are based, provide valuable insights for evaluating the fracture connectivity of a rock medium through which fluid can find stress-free pathways to flow, and the medium becomes a potential source of economic production. Fracture orientation, length and density play an indispensable role in the characterization of fracture networks and the estimation of connectivity in fractured media (Mauldon et al., 2001).

We investigated the connectivity of fracture networks using fracture geometric properties in this research. Specifically, fracture directions, the coefficient of dispersion, density, and fracture length were analyzed to determine the connectivity of fracture networks. The characterization of these fracture properties was demonstrated to evaluate the connectivity of the fracture system at an outcrop scale in 2D. The variation in fracture strike angle was measured by calculating the coefficient of dispersion (i.e., Fisher coefficient). The coefficient of dispersion depicts the distribution of fracture azimuth in the field area. Fracture sets observed on or near geological structures with a low Fisher coefficient value indicate a non-uniform orientation distribution, whereas a high value indicates a unidirectional fracture orientation pattern. A well-developed cluster is defined by the combinatorial relationship among strike data, FCA (a fraction of the fracture network’s connected area) and FCL (the distribution of trace length in terms of cluster size in the overall sampling area). Fracture length develops a negative exponential relationship with connectivity and has little effect directly on the connectivity of individual fractures, but greatly contributes to the interconnection of fracture networks within the fracture system. The determination of the size of clusters can have a substantial and influential effect on fracture connectivity in terms of the fractional connected area of the fractures in the sampled area at a given observational scale. The results of twenty-seven fracture trace maps (of which only three maps are presented in this study due to paper length constraints) show that a larger observational scale results in longer fracture trace length, implying a proportional correlation between trace length and observational scale and an inverse relationship with fracture density. The only scale-independent measure is fracture orientation. The fracture properties in a fracture system provide mutually predictable information about the connectivity of fracture networks at various scales.

Our findings, in this article, are consistent with Odling’s (1997) findings and validate the methodology used by Ghosh and Mitra (2009). It is also noted, as shown in the examples of the trace maps described above (Fig. 10, Fig. 11, Fig. 12, Fig. 13, and Fig. 14), that the number of fractures affects the density of fractures linearly, and that larger clusters are further connected by the backbone (longer) fractures, allowing fracture networks to be interconnected. Increased network interconnection makes the fracture system more productive in terms of easy fluid migration pathways.

This study does not include the characteristics of fracture aperture, the effects of lithological contrast on fracture density or the controlling factors of length, related to the deformation history of the structural domain.

7. Conclusion

The analysis in our case thus concludes that connectivity increases directly with the number of short-length fractures, the intensity of fracture orientation dispersion, fracture density, and the scale of observation. However, longer fractures do not directly affect the density of fractures but contribute to the interconnectedness of clusters and fracture networks. The field data are in good agreement with the trace maps illustrated in the figures described above.

The values of the coefficient of dispersion are estimated between 0 and 1. In our analysis, the orientation of all the fracture sets except Fracture set-2 is constrained to the value of dispersion coefficient 1, implying that the maximum population of fractures represents fracture sets concentrated in their strike angles and exhibit an orthogonal relationship with one another. The orthogonal geometry of the fractures produces a “T” intersection, whereas a conjugate fracture system increases the number of intersection nodes. This geometrical relationship further indicates that the fracture sets are systematic and can help to increase the likelihood of fracture interconnection. Fracture Set 2 has a uniform distribution with a low coefficient and serves as a “backbone” in fracture network connectivity. Fracture connectivity is essentially determined by the number of fracture sets, orientation dispersion, and fracture set density. Future research needs to incorporate these parameters to investigate how fluid transport and fractures interact in hydrocarbon reservoirs or aquifers. Using the Khairi Murat area as a case study and applying the findings and methodology from this study, it is possible to model the flow issues and transport characteristics of subsurface fractured reservoirs by learning more about the geometrical and hydraulic properties of fractures.