Multiscale fracture networks and their impact on hydroshearing response in the Canadian Shield (Kuujjuaq, Canada)

Abstract

Understanding the natural fracture network is essential for geothermal-related investigations. However, the geometrical attributes depend on the scale of observation. Therefore, a multiscale characterization of the fracture network is essential to ensure that heat and flow simulations are based on stochastically generated discrete fracture network models representative of the natural fracture system observed. The objective of this project was to understand the scale effect of fracture data on the performance of a potential enhanced geothermal system in crystalline rock in northern Canada. This was accomplished by collecting and characterizing fracture data from core, outcrops and satellite image, and then constructing a discrete fracture network model which was used to simulate the performance of the geothermal system. The numerical simulations suggest that fracture length and spacing have an important impact on its performance. Thermal short-circuiting can be easily achieved if the fracture network is modelled based solely on satellite image data, and hydroshearing may be less effective if the DFN is constructed based solely on outcrop data. The simulations suggest that combining the different datasets provides the best compromise between heat extraction, water losses, hydraulic impedance and thermal drawdown. Despite the uncertainties, the fracture data used highlights the importance of multiscale fracture analysis for heat-flow simulations of geothermal reservoirs.

Article Highlights

• A multiscale fracture network characterization using satellite imagery, outcrop and core data is presented.

• The scale effect of fracture properties on hydroshearing is explored.

• Fracture size and spacing play an important role in the development of engineered geothermal energy systems.

1 Introduction

In Canada, over 280 communities are disconnected from the provincial/territorial power grid and primarily relying on diesel power plants for electricity and on oil furnaces for space heating and domestic hot water (CER 2022). At the nowadays context where climate change and energy transition are buzzwords, many remote communities are exploring opportunities to reduce their reliance on fossil fuels, to improve energy security and to decrease the carbon footprint (CER 2022). One such opportunity includes developing renewables (CER 2022). In fact, a key component of the Pan-Canadian Governmental Framework (GCan 2022) is to reduce greenhouse gas emissions by helping remote communities transitioning from diesel-fired generation to cleaner energy sources. In this context, geothermal energy may be a viable option for certain communities.

According to the American Geosciences Institute (Neuendorf et al. 2005), a geothermal system is “Any regionally localized geological setting where naturally occurring portions of the earth’s internal heat flow are transported close enough to the earth’s surface by circulating steam or hot water to be readily harnessed for use”. Williams et al. (2011) broadened this definition, and for these authors “A geothermal system is any localized geological setting where portions of the Earth’s thermal energy may be extracted from a circulating fluid and transported to a point of use”. To include the concept of Enhanced or Engineered Geothermal Systems (EGS), Williams et al. (2011) definition can be extended to (Moeck 2014):

“A geothermal system is any localized geological setting where portions of the Earth’s thermal energy may be extracted from natural or artificially induced circulating fluids transported to a point of use. Enhanced Geothermal Systems are portions of the Earth crust where the ratio of flow rate and fluid temperature is naturally too low for economic use, and therefore the flow rate must be increased to a sufficient flow rate/temperature ratio by enhancing the natural permeability through technological solutions”.

An example of where EGS could be deployed is regions of crystalline bedrock, such as the vast Canadian Shield which covers most of northern Canada, where permeability is mainly supported by a network of natural fractures that may be open or sealed. In instances where the natural fracture network exhibits poor hydraulic connection, its permeability needs to be increased to permit fluid circulation through the rock mass at rates of commercial interest. Stimulation options include hydraulic, chemical, and/or thermal methods (e.g., Genter et al. 2010).

The potential of EGS to provide power and heat is enormous considering that at some depth, anywhere on earth, the rock is hot enough. However, there are a number of scientific and technical challenges to be overcome (e.g., Pollack et al. 2021), a key one being the ability to create a reservoir with sufficient permeability. This is greatly influenced by the geometry, connectivity, and frictional properties of the natural fracture network. This has been shown beginning with the first EGS field sites (e.g., Richards et al. 1994; Genter et al. 2010; Brown et al. 2012), and thus it is not surprising that detailed fracture characterizations have been completed at every geothermal site where EGS has been applied (e.g., Dezayes et al. 2010; Brown et al. 2012; Glaas et al. 2021; Reinecker et al. 2021).

Fracture properties can be measured at different scales, from the borehole (and core, ideally oriented), to outcrop, to images from planes and satellites. The scale of observation affects what fracture properties can be measured. For example, fracture strike and dip can be determined from borehole logging, oriented core samples, and outcrops, but cannot be measured from satellite images.

As well, for fracture properties that can be measured at all scales, the scale of observation has an influence on the value of these. For example, fracture density can be defined at a fine scale from borehole logging, oriented core, and outcrops, but only at a coarse scale from areal images. Another example is length: shorter fractures can be identified in outcrop logging, but larger fractures might only be detected (or at least their length quantified) using areal images.

Simulating the design and performance of an EGS involves the construction of a discrete fracture network model. The model is used to understand how the natural fracture system responds to stimulation treatments, and subsequently how the reservoir responds to the production and injection of fluids, and ultimately, how much thermal energy can be produced.

Since the scale of observation has an impact on the geometrical attributes of fractures, a multiscale fracture network characterization should be undertaken to avoid misrepresentation and mischaracterization of these attributes (e.g., Hardebol et al. 2015; Meixner et al. 2018; Espejel et al. 2020; Bossennec et al. 2021, 2022; Chabani et al. 2021a, b; Ceccato et al. 2022; Frey et al. 2022; Vega and Kovscek 2022). In fact, access to an exhaustive dataset, allowing mapping at different scales in both surface (e.g., outcrop mapping, drone imagery, airborne LiDAR) and subsurface (e.g., cores, borehole, seismic), seems ideal to provide a comprehensive setting to fully characterize fracture networks and additional structural properties (e.g., Espejel et al. 2020). Furthermore, beyond a comprehensive structural characterization of the naturally fractured geothermal reservoir, understanding the thermal, hydraulic and mechanical coupled processes during reservoir development and circulation is additionally important to forecasting EGS performance. This is done through numerical modeling, and several scholars have devoted themselves to understanding the influencing parameters that affect production performance (e.g., Chen et al., 2019; Ling et al., 2019; Zhang et al., 2019; Liu et al., 2021).

As aforementioned, ideally, fracture data for an EGS study site should be collected at all scales and incorporated into the DFN to create the most realistic depiction of the natural fracture network. However, this is not always feasible in remote northern communities. Many of these communities do not have boreholes, or cores recovered from wells, and in many cases the quality of the satellite images is poor for a good analysis. Moreover, not all communities are covered by LiDAR. Fortunately, Kuujjuaq, the study site has a relatively good dataset to study the scale effect on fracture properties. Kuujjuaq is an Indigenous community in Nunavik (northern Quebec) that relies entirely on diesel for its heat and power.

The objective of this paper was to analyze fracture properties at the core, outcrop, and satellite image scale in the vicinity of Kuujjuaq. Then, for each scale, construct a DFN and use it to simulate the response of the fracture network on hydroshearing and the performance of the stimulated reservoir. The numerical simulations permitted to conduct a qualitative analysis of how the scale effect could influence the performance of an EGS. These highlight how length and spacing could play a key role in the reservoir development, shedding light on the need for multiscale fracture data analysis to be used in simulating EGS.

The results of this research work can be used to guide further studies in remote northern communities. This work illustrates that simple and affordable methods such as outcrop and satellite imagery analysis can help to better understand how an EGS could perform in the Canadian Shield crystalline basement.

2 Geology

2.1 Local geological setting

Kuujjuaq is located on the Southeastern Churchill tectonic province of the Canadian Shield. Two lithologies are dominant in the study area – paragneiss and diorite (Fig. 1). Smaller outcrops of gabbro, tonalite and granite are also observed (Fig. 1). The paragneiss belongs to the False Suite (Archean), the diorite is part of the Kaslac Complex (Paleoproterozoic), the gabbro is associated with the Ralleau Suite (Archean to Paleoproterozoic), the tonalite belongs to the Aveneau Suite (Paleoproterozoic) and the granite to the Dancelou Suite (Paleoproterozoic).

Fig. 1

Geology of Kuujjuaq (adapted from MERN 2019). Black square on the geological map indicates the location of satellite image for structural lineament interpretation. LP: Lac Pingiajjulik fault, LG: Lac Gabriel fault, F: Forum well

Two main regional faults cross the community. One of these structures is the Lac Pingiajjulik fault (Fig. 2) and the other is an interpreted fault from analysis of aeromagnetic images by previous researchers. According to Poirier (1989), the Lac Pingiajjulik fault represents a thick zone of highly recrystallized mylonite that separates two lithological blocks of different ages However, aeromagnetic images suggest a different interpretation. The exact location of this fault is still debatable, and that discussion goes beyond the scope of this work. This structure is masked in several places by Paleoproterozoic intrusive rocks of the Aveneau and Dancelou suites, suggesting that the fault predates the emplacement of these intrusive units.

Fig. 2

Mylonitic foliation suggesting the presence of the Lac Pingiajjulik fault in horizontal (left-hand side) and vertical (right-hand side) surfaces

The Lac Pingiajjulik fault is associated with a NW-verging thrust phase, followed by a late dextral movement (Perreault and Hynes 1990; Clark and Wares 2004). Rocks located to the east of the Lac Pingiajjulik fault are metamorphosed to the upper amphibolite facies and granulite facies. Rocks located to the west of this structure are metamorphosed to the upper amphibolite facies. Strong recrystallization of mylonite indicates that metamorphic conditions were sufficient to allow static recrystallization after the main episode of deformation (Poirier et al. 1990). The Lac Pingiajjulik fault corresponds to a more or less sharp lineament that follows the contour of lithodemic units on the regional aeromagnetic map (MERN 2020).

2.2 Natural seismicity

The Canadian National Seismic Network has a station in Kuujjuaq since 1998. Only three earthquake records occurring more than 10 km away from the community of Kuujjuaq were found in the Canadian database (GCan 2021):

• (1) Event ID 20111118.1859002 had a magnitude of 2.2, a depth of 18 km and occurred 10 km N from Kuujjuaq;

• (2) Event ID 20160313.0434004 had a magnitude of 3.4, a depth of 18 km and occurred 53 km N from Kuujjuaq;

• (3) Event ID 20160313.0435001 had a magnitude of 3.0, a depth of 18 km and occurred 36 km NW from Kuujjuaq.

This suggests that the study area is placed in a aseismic zone of the Canadian Shield craton.

3 Description of the datasets and data analysis

Structural lineaments were analyzed by interpretation of Google Earth satellite images and fracture information (i.e., azimuth, dip, length, spacing, frequency, aperture and connectivity) was collected in outcrops and well cores.

A structural lineament was defined based on O’Leary et al.'s (1976) definition of the word “lineament”: “A lineament is a mappable, simple or composite linear feature of a surface, whose parts are aligned in a rectilinear or slightly curvilinear relationship, and which differs distinctly from the patterns of adjacent features and presumably reflects a subsurface phenomenon.”

3.1 Satellite imagery

The satellite imagery for interpretation of the structural lineaments was retrieved from Google Earth (Fig. 3). The outcrop analyzed is located on the east side of the Koksoak River and it is made mostly of gneissic rocks of the Ungava Complex (see Fig. 1 for geographical location). This image was chosen since it was the one with the best quality for analysis. The images near the community have a lower quality for interpretation of the structures.

Fig. 3

Google Earth satellite image and rose diagram produced from structural lineament interpretation

The satellite image was interpreted in ArcGIS version 10.3 using the Editor—straight element function to manually depict the structural lineaments. Then, trace azimuth and length were evaluated as explained in the ESRI technical support webpage—How to: Batch calculate line direction using the Field Calculator (ESRI 2021). The information gathered about trace orientation was then used in Grapher (Golden Software 2022) to produce the rose diagram.

A stochastic DFN model is commonly built on the statistical information provided by the fracture data analysis. A general stochastic DFN approach assumes fractures to be straight lines or planar discs and treats the other geometrical properties as independent random variables obeying certain probability distributions (Lei et al. 2017). For this reason, probability distribution functions (PDFs) were applied to characterize the trace azimuth and length (Appendix 2) and statistical tests such as the root mean squared error (RMSE) and the sum of squared errors (SSE) were used to evaluate which PDF best fit the data sampled (Appendix 3) as suggested by Zeeb et al. (2013b).

A total of 226 structural lineaments were identified, suggesting an aerial frequency of 13.29 km⁻². Interpretation of the structural lineaments suggests a main lineament orientation E-W and NE-SW (Fig. 3).

3.2 Outcrops

Six sampling areas were selected to carry out fracture analysis (A1 to A6; see Fig. 1 for geographical location). Sampling area 1 (A1) is located near the Coop store in Kuujjuaq (Fig. 4). In this area, a total of 6 outcrops were analyzed (Appendix 1). Sampling area 2 (A2) is located outside of the community in an area informally called Trois lacs (Fig. 4). In this area, 5 outcrops were investigated (Appendix 1). Sampling area 3 (A3) is located within the community, near a new treatment center and Koksoak River (Fig. 4). In this area, fracture data was collected in 8 outcrops (Appendix 1). Sampling area 4 (A4) is located near the Hydro-Quebec diesel power plant (Fig. 4). In this area, 3 outcrops were analyzed (Appendix 1). Sampling area 5 (A5) is located near the graveyard (Fig. 4). In this area, fracture data was sampled in 3 outcrops (Appendix 1). Sampling area 6 (A6) is located near the dock (Fig. 4). This is the only outcrop with a vertical surface.

Fig. 4

Representative photos of the different fracture sampling areas analyzed and rose diagrams and stereonet produced from fracture interpretation. Red line in the stereonet indicates the outcrop orientation

The outcrops in sampling areas A1 and A2 are made of paragneiss rocks from the False Suite and the outcrops in A3 and A4 are made of diorite rocks from the Kaslac Complex. The outcrops in sampling area A5 are near the contact between the paragneiss rocks of False Suite and tonalite rocks of Aveneau Suite and fractures from both rock types were investigated. Sampling area A6 corresponds to the contact between gabbro rocks from the Ralleau Suite and tonalite rocks from the Aveneau Suite, and fractures from both lithologies were investigated.

The linear scanline sampling method together with a transit compass corrected for the magnetic declination (22° in Kuujjuaq) and an electronic caliper were used to collect fracture geometrical attributes (i.e., azimuth, dip, size or length and aperture) and the position of fractures along the scanline. A total of 1 to 4 scanlines per outcrop were used in an attempt to fully characterize the outcrop surface. The number of scanlines per outcrop was chosen depending on the outcrop quality and homogeneity/heterogeneity of fractures observed.

Grapher (Golden Software 2022) was used to build the rose diagrams considering the azimuths of the fractures observed and Stereonet version 11.3.7 (Allmendinger 2020) to produce stereonet canvas. PDFs were applied to characterize fracture properties (Appendix 2). Then, statistical tests such as the RMSE and the SSE were used to evaluate which PDF best fit the data sampled (Appendix 3) as suggested by Zeeb et al. (2013b).

3.3 Core

In September 2021, a 240 m deep well was drilled near Kuujjuaq’s Forum building (Fig. 1) to support investigations into the potential to use a ground-source heat exchange system to provide heat for the building (Geotherma solutions Inc, 2022). This well also helped to further advance research into this community deep geothermal potential, and geothermal energy development, including EGS. The entire length of the well was cored.

The core was logged for lithology and fractures (Fig. 5). The main lithology encountered was paragneiss rocks from the False Suite, interlayered with diorite and gabbro rocks from Kaslac Complex and Ralleau Suite and granitoid rocks from the Kaslac Complex and Aveneau and Dancelou suites.

Fig. 5

Graphical well log and representative photos of the lithologies encountered and natural fractures identified. NF: possible natural fractures displayed in stick plot and absolute frequency histogram for 3 m intervals

The linear scanline sampling method was used to locate the fractures along the well core and obtain information about fracture frequency and spacing. The tape was placed parallel to the core length and thus vertical structures are underrepresented since these are parallel to the scanline. PDFs were applied to characterize the data gathered (Appendix 2) and the same statistical tests mentioned above were used to evaluate which PDF best fit the data sampled (Appendix 3) as suggested by Zeeb et al. (2013b).

A total of 187 possible natural fractures were identified, suggesting a linear fracture frequency of 0.80 m^− 1 and a mean apparent fracture spacing of 1.11 m (Fig. 5). The interpretation took into consideration the shape of the fractures and the mineralization. The core was not oriented and borehole images were not taken (e.g., borehole televiewer) and therefore information about fracture orientation and dip was not retrieved. Further information about geological core logging is given by Miranda and Raymond (2022). The interpretation of the well core suggests neither a relationship between lithology and fracture frequency, nor between fracture frequency and depth (Fig. 5).

Histogram is the most common tool to present fracture data results. However, binning the data in fixed intervals may degrade the information. For that reason, Sanderson and Peacock (2019) recommend the use of stick plots, for example, to preserve more information about individual fractures. Stick plots provide a good visual representation of the fracture distribution along the line, while this information is lost by using histograms (Fig. 5).

4 Discrete fracture network model

4.1 Influence of lithology

The first fracture data analysis carried out had the goal of understanding if lithology influences the geometrical properties of the fractures observed at the outcrop scale and consequently hydroshearing.

Sampling areas A1, A2, A3 and A4 were analyzed with this objective in mind. As mentioned in Sect. 3.2., sampling areas A1 and A2 are located on paragneiss rocks while the outcrops in A3 and A4 are made of diorite rocks.

The results suggest that the paragneiss rocks have three main fracture sets given by a Mixture of von Mises (MvM) distributions with 18 elements (Fig. 6): F1—N81-90°E (E-W), F2—N151-160°E (NW–SE) and F3—N1-10°E (N-S). Three main fracture families were also identified for the diorite rocks which are (Fig. 6): F1—N91-100°E (E-W), F2—N121-130°E (NW–SE) and F3—N151-160°E and N161-170°E (NNW-SSE). These results seem to suggest that, despite the different distribution of values, both lithologies have roughly the same fracture sets.

Fig. 6

Influence of lithology on the geometrical attributes of fractures

Fracture size is also similar for both paragneiss and diorite (Fig. 6). These lithologies revealed a median fracture size of 1.4 m and 2 m, respectively for paragneiss and diorite. However, the diorite rocks have a larger dispersion of fracture size than paragneiss rocks.

Similarly, minimal differences were also found for fracture spacing (Fig. 6), with the paragneiss revealing a mean fracture spacing of 0.18 m and the diorite 0.36 m, with the latter showing larger data dispersion.

Aperture data also reveals minimal differences (Fig. 6), with both paragneiss and diorite revealing the same median aperture of 1.5 mm and similar data dispersion.

4.2 Influence of scale effect

The structural lineaments and fracture data collected from the satellite image, outcrops and well core were analyzed separately and combined to, first, understand how the scale of observation influences the geometrical attributes of fractures, second, evaluate which probability distribution function (PDF) best fits the data collected, and third, investigate the influence of scale effect on hydroshearing.

A MvM distributions with 18 elements was applied to the outcrop data and suggests a unimodal distribution, with one single main fracture set oriented in the E-W direction (N81-90ºE; Fig. 7). Analysis of the satellite image, on the contrary, suggests a multimodal distribution with five main fracture sets (Fig. 7): F1—N81-90°E (E–W), F2—N41-50°E (NE–SW), F3—N101-110°E (NW–SE), F4—N21-30°E (NNE–SSW) and F5—N151-160°E (NNW-SSE). Combining both datasets results also in a multimodal distribution, but the set E-W is the dominant one (Fig. 8).

Fig. 7

Statistical analysis of the geometrical attributes

Fig. 8

Influence of lithology on the shear displacement

Analysis of the size suggests that 50% of the interpreted structures in the satellite image are above 300 m, while 50% of outcrop fractures are below 2 m (Fig. 7). Combining both datasets suggest that the DFN would be characterized by 65% of structures below 50 m and 35% of structures between 50 to 1,437 m. This highlights the greater number of small fractures observed compared to larger structures, agreeing with the fact that a power-law distribution provides the best fit to the structures size. The RMSE for the power-law distribution was evaluated to be 0.25, while for the exponential and lognormal distribution was 0.29 and 0.35, respectively. The SSE for the power-law distribution was estimated to be 40, while for the exponential and lognormal distribution was 52 and 80, respectively.

As the power-law distribution best describes the fracture sizes, it can be assumed that the network follows a fractal size distribution. Plotting the data collected in outcrops and observed in the satellite image in a log–log plot suggests a multifractal size distribution and application of the chord method suggests a lower cutoff fracture size of 200 m when both datasets are combined. The power-law exponent for fractures between 200 to 1437 m is −2.0. The power-law exponent for the outcrop data is −1.8 and the lower cutoff length is 2 m. The power-law exponent for the satellite image data is −2.0 and the lower cutoff length is 200 m. The chord method is based on the log–log plot of fracture size against cumulative frequency. A line is drawn through the data point with the shortest length and the data point with the longest length. Then, the data point with the longest perpendicular distance from this line is used as the lower cutoff (Zeeb et al. 2013a).

The power-law functions to characterize our dataset are then:

$${\rm Outcrop\;data}:f\left( l \right) = 1.407l^{ - 1.792}$$

(1)

$${\rm Satellite\;image\;data}:f\left( l \right) = 39982l^{ - 1.982}$$

(2)

$${\rm Both\;databases}:f\left( l \right) = 14230l^{ - 1.982}$$

(3)

These power-law exponents closer to 2 suggest that the fracture size has a fractal character. Furthermore, the power-law exponent of the outcrop data suggests that based on that dataset the fracture network comprises short size elements (Bossennec et al. 2021).

Spacing data interpreted suggests that the best fit is given by an exponential distribution, which indicates a uniform spatial distribution for the structures (Fig. 7). The RMSE for exponential distribution is 0.06 and the SSE is 2. The power-law distribution, on the other hand, revealed RMSE and SSE of 0.56 and 159. Furthermore, analysis of the data indicates that the core has more widely spaced fractures than the outcrops. Average fracture spacing obtained from the core data is 1.11 m while outcrop average fracture spacing is 0.55 m. This agrees with the linear fracture frequency estimated. Core data revealed a value of 0.80 m⁻¹ while outcrop data suggest an average value of 1.81 m⁻¹. In other words, while core data indicates less than a fracture per meter of core, outcrop data implies almost 2 fractures per meter of outcrop.

Fracture aperture measured in the field with an electronic caliper indicates a mean value of 6.5 mm, with a minimum value of 0.1 mm and a maximum value of 165 mm (Fig. 7). These values indicate that the apertures range from tight to wide, but with different probability of occurrence (Table 2).

The data collected also suggests that a power-law distribution has the best fit, with an RMSE of 0.17 and an SSE of 9. A lognormal distribution reveals an RMSE of 0.29 and a SSE of 27.

The power-law function to characterize our dataset is then:

$${\rm Outcrop\;data}:f\left( w \right) = 0.3322w^{ - 0.3668}$$

(4)

It is important to highlight, however, that these apertures are greatly influenced by weathering of the fracture walls and may not correspond to the in situ aperture value. For this reason, these aperture values measured in the field were neglected for the numerical simulations undertaken.

5 Description of the numerical model

As described in the Introduction section, numerical simulations are used to investigate the development and performance of EGS.

The construction of a geological model that represents the natural fracture system (and includes other rock properties as appropriate to the simulation goals) is the first step.

FRACSIM3D is a fracture generator and thermo-hydro and hydro-mechanical coupled model (Willis-Richards et al. 1996; Jing et al. 2000; Miranda et al. 2021). We used FRACSIM3D to construct the geological model. FRACSIM3D can generate fracture networks by taking the fracture attribute data, and the PDFs for each fracture attribute, and generates a natural fracture network. This numerical simulator was also used to simulate the hypothetical EGS.

A detailed description of the numerical model (i.e., fracture generation process, fracture constitutive laws implemented, and hydroshearing, circulation and heat extraction processes) is given by Willis-Richards et al. (1996) and Jing et al. (2000).

The validity of using FRACSIM3D is supported by the work of Jing et al. (2000) who used FRACSIM3D to simulate the short- and long-term thermal extraction of an EGS at Hijiori (Japan) as the numerical model results were comparable with the field tests.

For our study site, Miranda et al. (2021) used FRACSIM3D to simulate the fracture sets that were optimally oriented to slip, and the results were similar to a method based on Mohr–Coulomb analysis carried out by Miranda et al. (2023). Furthermore, Miranda et al. (2021) used FRACSIM3D to investigate the fluid pressure necessary to induce slippage along fractures, and the results were similar to the estimates carried out by Miranda et al. (2023) using the critical fluid pressure estimates method.

5.1 Conceptual model

A cubic model volume of 8 km³ (with 2 km of edge length) discretized into a grid of 200 × 200 × 200 cells was used to carry out the numerical simulations. A grid dependency study was undertaken and revealed a difference of about 5% between a 100-cell grid and a 300-cell grid. The fracture network is incorporated into the model. The system was designed as a doublet, with one injector and one producer, both vertical. A spacing of 600 m was assumed between the wells and an open hole length of 600 m was considered. The wells were positioned oriented parallel to the maximum horizontal principal stress. The well bore radius was defined as 0.11 m. A stimulated rock volume of 0.8 km³ was defined for both wells.

5.2 Fracture constitutive laws

The following equations implemented in FRACSIM3D are intended to be practical approximations to the shear dilation behavior of a model containing thousands of fractures.

5.2.1 Initial fracture aperture

Fractures in FRACSIM3D are assumed to be frictional contacts without any infilling minerals causing resistance to sliding or opening (Jing et al. 2000). The initial fracture apertures, at zero effective stress, are assumed in FRACSIM3D to be proportional to the fracture radii (Willis-Richards et al. 1996; Jing et al. 2000):

$$a_{0} = \beta r_{{\text{i}}}$$

(5)

where a₀ (m) is the undisturbed fracture aperture and β (-) is a constant of proportionality chosen such that the macroscopic permeability of the undisturbed network approximately matches the in situ measured permeability. This constant of proportionality is uniformly applicable to all the fractures in the network and is calculated by (Willis-Richards et al. 1996):

$$\beta = \left( {\frac{{12Sk\mathop \sum \nolimits_{{\text{all i}}} r_{{\text{i}}} }}{{\mathop \sum \nolimits_{{\text{all i}}} \frac{{r_{{\text{i}}}^{4} }}{{\left( {1 + 9\frac{{\left( {\sigma_{{\text{n}}} - P_{{\text{p}}} } \right)}}{{\sigma_{{{\text{n}},{\text{ref}}}} }}} \right)^{3} }}}}} \right)^{\frac{1}{3}}$$

(6)

where S (m) is the mean fracture spacing measured in the field, k (m²) is the mean in situ permeability from field measurements, σ_n (Pa) is the rock stress normal to the fracture surface, P_p (Pa) is the fluid pressure in the fracture and σ_n,ref (Pa) is the effective normal stress applied to cause 90% reduction in the compliant aperture.

5.2.2 Fracture shear dilation behavior

Fractures will begin to slip when the shear stress acting parallel to the fracture plane exceeds its sliding resistance. The shear stability of a rough fracture is expressed in FRACSIM3D through the linear Mohr–Coulomb criterion (Jing et al. 2000; Miranda et al. 2021):

$$\tau > F_{{{\text{asperity}}}} \left( {\sigma_{{\text{n}}} - P_{{\text{p}}} } \right)tan\left( {\phi_{{{\text{basic}}}} + \phi_{{{\text{dilation}}}}^{{{\text{effective}}}} } \right)$$

(7)

where τ (Pa) is the shear stress, F_asperity (-) is the asperity strength factor, ϕ_basic (º) is the basic friction angle of the rock material and ${\phi }_{\mathrm{dilation}}^{\mathrm{effective}}$ (º) is the effective shear dilation angle. This property is related to the shear dilation angle at zero normal stress, and is expressed approximately in FRACSIM3D as (Jing et al. 2000):

$$\phi_{{{\text{dilation}}}}^{{{\text{effective}}}} = \frac{{\phi_{{{\text{dilation}}}} }}{{1 + 9\frac{{\left( {\sigma_{{\text{n}}} - P_{{\text{p}}} } \right)}}{{\sigma_{{{\text{n}},{\text{ref}}}} }}}}$$

(8)

Shear displacement is expressed in FRACSIM3D as (Jing et al. 2000):

$$U = \frac{{\tau - \left( {\sigma_{{\text{n}}} - P_{{\text{p}}} } \right)tan\left( {\phi_{{{\text{basic}}}} + \phi_{{{\text{dilation}}}}^{{{\text{effective}}}} } \right)}}{{K_{{\text{s}}} }}$$

(9)

where U (m) is the shear displacement and K_s (Pa) is the fracture shear stiffness, which is approximated in FRACSIM3D as (Jing et al. 2000):

$$K_{{\text{s}}} = \gamma \frac{G}{r}$$

(10)

where γ (−) is a geometric parameter ($\gamma =\frac{7\pi }{24}$ for circular fractures; Eshelby 1957) and G (Pa) is the shear modulus of the rock material.

The aperture of a sheared fracture with unabrabed asperities in contact is then approximated in FRACSIM3D as (Jing et al. 2000):

$$a = \frac{{a_{0} + U \times tan\left( {\phi_{{{\text{dilation}}}}^{{{\text{effective}}}} } \right)}}{{1 + 9\frac{{\left( {\sigma_{{\text{n}}} - P_{{\text{p}}} } \right)}}{{\sigma_{{{\text{n}},{\text{ref}}}} }}}}$$

(11)

Rock mechanics experiments have suggested that shear dilation angle increases over the first few millimeters of displacement as asperities of increasing wavelengths become out of phase, reaching a maximum and then decreasing to lower values after few millimeters of displacement since the longer wavelength asperities have less relative amplitude than those with the shorter wavelength. Thus, a shear displacement–dilation relationship, in which the shear dilation angle is a reducing function of displacement rather than constant, was formulated and implemented. The relationship is as follows:

• A linear increase from a starting shear dilation angle value to a maximum value over a certain small shear displacement distance;

• An exponential decay with displacement thereafter to a low constant value at greater displacement at a user-specified rate.

5.3 Model assumptions

The reservoir simulated was assumed homogeneous and isotropic in terms of thermal, hydraulic and mechanical properties. The mechanical properties assigned to each fracture (i.e., fracture basic friction angle, shear dilation angle, and closure properties) in the model were assumed to be the same to simplify the modeling process. Assumptions for the properties of the medium are taken from Miranda et al. (2021; Table 1). The frictional properties of fractures used in this study are working hypotheses. The simulations assumed that the center of the reservoir is placed at 4 km depth.

Table 1 Medium properties assumed as input parameters for simulations (Miranda et al., 2021)

5.4 Boundary and initial conditions

The initial reservoir temperature at the center of the model is defined by the user, and in this case the value 88 ºC was used. The upper and boundary temperature conditions are calculated based on this temperature value and the geothermal gradient. No heat transfer occurs beyond the boundaries of the stimulated rock volume.

The initial reservoir pressure at the center of the model is assumed hydrostatic, and the value used for the in situ fluid pressure at model center is 43 MPa. The model is assumed to be filled with water. The upper and boundary pressure conditions are calculated based on this pressure value and the pressure gradient.

The initial reservoir permeability is defined by the user, and a value of 10^–17 m² was used in this case. The initial fracture aperture is calculated numerically considering this initial in situ permeability and magnitude of the principal stresses imposed in the model, and the numerically generated stochastic fracture network.

A larger permeability is assigned to the far field (i.e., beyond the boundaries of the model) to permit flow out of the model through fractures that do not cross the production well. This condition also permits flow to be brought into the model by fractures whose length exceeds the model boundaries. The permeability of the far field was increased by a factor of 10 in these simulations.

As aforementioned, a stimulated rock volume surrounding the injection (and production) well(s) is defined by the user. In this case, we assumed that both wells were stimulated, and the stimulated rock volume was 0.8 km³ in both wells. This stimulated rock region constrains the volume of rock containing fractures that will be disturbed by changes imposed into their initial state of stress. Beyond this volume boundaries, the pressure drops rapidly towards initial conditions. Similarly, within this stimulated volume, the pressure drops by 25% from the center of the stimulated region, where the well openhole is, towards the edge of the pressurized region. The user defines the stimulation fluid pressure at the wells. A value varying between 47 and 51 MPa was used in the simulations. These stimulation fluid pressures were previously calculated as necessary to reactivate the optimally oriented fractures through slip tendency analysis (Miranda et al. 2023).

This change into the fractures initial state of stress leads to the aperture enhancement of the fractures that are well oriented considering the direction of the maximum horizontal principal stress imposed by the user (N215ºE, in this case; Miranda et al. 2023). If a fracture falls inside of the stimulated region, it will be reactivated. If the fracture falls inside and outside of the stimulated region, this fracture will be reactivated. If the fracture falls outside of the stimulated region, it will not be reactivated.

During the circulation and heat extraction phases, a pressure varying between 14 and 30 MPa was imposed to the injection well and a pressure of −2 to −20 MPa was imposed to the production well for the stimulation scenarios. A fluid temperature of 30 °C was imposed in the injection well. Heat is extracted from the reservoir through the reactivated fractures. The simulations were run for a period of 20 years.

6 Numerical simulation results

6.1 Influence of lithology

The numerical simulations suggest that lithology has a minimal influence on the shear displacement (Fig. 8). The median shear displacement for the paragneiss is 0.99 mm and 0.98 mm for the diorite. The diorite shows greater dispersion of the shear displacement than the paragneiss.

Considering these observations, we concluded that lithology has a negligible effect on the geometrical properties of fractures and on the hydroshearing treatment. For this reason, these two lithologies were grouped into one single hydrostratigrahic unit for the following data analysis and modeling.

6.2 Influence of scale effect

The scale effect has a significant impact on the hydroshearing treatment (Fig. 9). The structural lineaments depicted from the satellite image reveal greater shear displacement than the fractures mapped on the outcrops. The median shear displacement is 5 mm for the satellite image data and 1 mm for the outcrop data. The satellite image data shows greater dispersion than the outcrop data.

Fig. 9

Influence of scale effect on the shear displacement

One of the causes for these differences may be related to the fracture shear stiffness. This parameter is approximated in FRACSIM3D as the ratio of the rock shear modulus divided by the fracture radius and multiplied by a geometric parameter. For the same shear modulus, a fracture with 0.08 m has a shear stiffness of about 333 GPa m⁻¹ while a 1437 m fracture has 0.02 GPa m⁻¹ of shear stiffness. The lowest shear stiffness values are thus found for the fracture network made of large fractures (satellite image) while the greatest values for the fracture network made of small fractures (outcrop).

These results suggest that hydroshearing treatment may be more effective in a network made of larger fractures than in a network made of smaller fractures. Larger fractures may have smaller shear stiffness and therefore shear displacement may be greater in these structures. This ultimately may have an impact on the sheared aperture, with larger fractures having larger void space than smaller fractures (Table 2).

Table 2 Aperture classification (ISRM 1978) and probability of occurrence

6.3 Impact on hypothetical EGS development and performance

As concluded in the above section, the length of the fractures could have an important impact on EGS development. A DFN built based exclusively on satellite image data could lead to an early thermal short-circuiting of the reservoir since only a few fractures may be fluid flow conduits, contributing to fluid circulation within the reservoir where only a few structures will act as preferential fluid flow pathways. This hypothesis was tested in this section by simulating nine scenarios where fracture size and frequency were changed but fracture orientation was kept constant (Table 3). The orientation of fractures observed from outcrop mapping was preferred since not all structural lineaments depicted in the satellite image may correspond to fractures.

Table 3 Numerical simulation scenarios

For the structural lineaments extracted from the satellite image, a slip patch radius of 150 m was assumed as a working hypothesis and to prevent unrealistic flow rates. In scenario 2C, the pressure at which the wells are held during circulation was modified to 30 MPa (injector) and -20 MPa (producer) and in scenario 2D the stimulation pressure in both wells was increased to 51 MPa. Scenario 3A corresponds to the untruncated fracture size while scenarios 3B and 3C were built considering the fracture size within the 30–70% interval. This fracture size range was chosen since fractures commonly consist of many interacting and linked segments, and what appears to be a single structural lineament mapped from satellite image may in fact be several structures (Peacock et al. 2019).

The simulations suggest that DFNs built exclusively based on outcrop data (Scenario 2A–2C) could lead to a greater reservoir impedance and lower heat extraction compared to a DFN built exclusively based on satellite image data (Scenario 1A–1B) or a mix of outcrop and satellite image data (3A–3C; Fig. 10). However, a DFN built exclusively based on the satellite image data may lead to an early thermal short-circuiting of the reservoir, as only a handful of fractures will be conducting flow (Scenario 1A; Fig. 10). Nevertheless, this behavior changes if the density of the network is increased (Scenario 1B; Fig. 10), or if outcrop data is introduced in the network (Scenario 3A–3C; Fig. 10). An interesting point illustrated by the simulations undertaken is that a DFN built mixing satellite image and outcrop data may need greater fracture density to achieve the same flow rates as a DFN constructed based on the satellite data alone. Water losses are also decreased in this mixed DFN with greater fracture density. Furthermore, a DFN constructed based solely on outcrop data may be more productive with greater fluid stimulation and circulation pressures.

Fig. 10

Top view slice at the center of the model illustrating the influence of fracture size, fracture frequency and fluid stimulation and circulation pressure on the development of an engineered geothermal system. White dot: injection well, red dot: production well

Additionally, the capacity of the EGS simulated to fulfil the heating demand of Kuujjuaq was investigated. The annual average heating needs of a typical residential dwelling in Kuujjuaq have been evaluated to be 22–71 MWh (Yan et al. 2019; Gunawan et al. 2020). This community has 973 dwellings, implying an annual average heating demand of 21 to 69 GWh. Furthermore, an annual population growth of 3% has been noticed through Census. Thus, assuming the same growth rate in heating demand, in 20 years, the annual average heating demand of this community may be 38 to 121 GWh. Keeping these values in mind, the simulation results indicate that Scenario 1A, 1B, 2C, 2D, 3A and 3C are able to fulfil the demand over the lifespan of the system considering the lower bound of demand (Table 3). However, these conclusions change considering the upper bound of heating demand. In this case, only Scenario 1B has sufficient thermal extraction to fulfil the heating demand (Table 3). However, the best compromise between heat extraction, water losses, reservoir hydraulic impedance and thermal drawdown is given by Scenario 3C, where a mix between small and large fractures spaced by 5 m was considered.

However, it is important to highlight that the results of these numerical simulations in terms of flow rate and thermal energy extraction are purely hypothetical and may differ once an EGS is developed. Field tests are necessary to validate the modeling results. The numerical simulations were carried out with the goal of illustrating qualitatively, rather than quantitatively, how different methods to characterize fracture length and fracture frequency may impact significantly EGS development and performance.

7 Discussion

Multiscale fracture network characterization has been investigated by many scholars (e.g., Hardebol et al. 2015; Meixner et al. 2018; Espejel et al. 2020; Bossennec et al. 2021, 2022; Chabani et al. 2021a, b; Ceccato et al. 2022; Frey et al. 2022; Vega and Kovscek 2022). All these studies highlight differences in fracture orientation, fracture length, and fracture spacing as a consequence of the scale effect. Since geometrical properties of fractures are scale-dependent, we decided to carry out the same investigation to quantify these differences using the crystalline basement in Kuujjuaq as a case-study. This multiscale characterization of the geometrical properties was then used to understand and gain insights on how the scale effect can influence mechanical simulations of hydroshearing treatment, and consequently, the potential development and the performance of a hydraulically stimulated geothermal reservoir.

The analysis of multiple data sources, namely satellite imagery, outcrops and core, suggests that fracture orientation, fracture length, and fracture spacing are indeed dependent on the scale of observation. Structural lineaments analyzed from the satellite image revealed a more complex distribution of orientations than the outcrop field analysis. One of the reasons may be due to the fact that not all structural lineaments observed from the satellite image may correspond to open-mode fractures. For that reason, the orientations retrieved from the satellite imagery analysis were considered highly uncertain and neglected in the numerical simulations carried out. Using outcrop data to model subsurface data is also highly uncertain, and data from core and well logs should be prioritized. Unfortunately, the well available for this study was drilled with other goals than geological data acquisition, and logging data was not acquired. Despite the availability of core samples, the core was not oriented, and the structures analyzed were considered only probable natural fractures since well logs would be needed to confirm if they are actually natural or rather induced fractures. In this work, we assumed the outcrop data to be the most accurate, but we are aware that we are introducing important bias and uncertainty to our results. Nevertheless, this study can inform about questions that need to be resolved and can trigger interest for further developments.

Furthermore, structural lineaments from the satellite image revealed larger fracture length than observed at the outcrop scale. This is expected since fractures are often not continuous segments, but rather a combination of interacting and linked segments as pointed out by Peacock et al. (2019). A large structural lineament observed at the satellite image scale can indeed correspond to several smaller fractures at the outcrop scale. This is an important data resolution bias and introduces important uncertainty to conduct further simulations with subsurface numerical models. As the numerical simulations illustrate, a DFN built based on satellite image data, where a majority of fractures are longer than 200 m, could lead to an early thermal short-circuiting of the geothermal reservoir. This can have important economic impacts for a geothermal system. However, if the DFN is built based on outcrop data, where most of the fractures have a length of 5 m or less, this could have important technological issues in developing the geothermal reservoir by hydraulic stimulation. In fact, the effectiveness of the hydraulic stimulation treatment may be greatly diminished with such a dataset, according to the numerical simulations carried out. Of course, field tests are needed to corroborate or refute the modeling exercise undertaken.

One option to deal with this data resolution bias according to Peacock et al. (2019) is to define a resolution limit, and sampling only fractures of a certain length. According to these authors, defining this resolution limit is helpful to upscale the geometry of the fracture network to the scales of interest for reservoir models. In this paper, we decided not to sample only fractures of a certain length since our objective was to understand how fracture lengths measured at different scales would respond to hydroshearing.

It is also important to discuss the fact that the satellite image analyzed corresponds to an outcrop made of gneiss of the Ungava Complex, while the field analysis was undertaken on outcrops of paragneiss of the False Suite unit, diorite of the Kaslac Complex and on outcrops on the contact between tonalite of the Aveneau Suite and paragneiss of the False Suite and between the tonalite and the gabbro rocks of the Ralleau Suite. Furthermore, the core is made of the paragneiss, diorite and tonalite units. The fact that the analysis groups different rock types can additionally create biases and uncertainty related to geology, as explained by Peacock et al. (2019). Despite the fact that the outcrop data suggests an independence of geometrical properties of fractures and geology, there is always a certain degree of geology bias involved. Ideally, the satellite image should have been made over paragneiss and diorite to decrease this bias. Unfortunately, the quality of the images available in Kuujjuaq limits this ideal situation, and we opted for the gneiss outcrop to illustrate the goal of this paper, which was to understand and gain insights on how the scale effect can influence the hydroshearing treatment and the potential development of an EGS in this remote northern community. This problem can be commented on to remote regions where subsurface data is sparce and expensive to collect.

Methods can also impose biases and uncertainty as explained in Peacock et al. (2019). Despite different methods being used in this work, namely analysis of a satellite image, analysis of core data, and scanline sampling methods on outcrops surface, the dataset used in this work was useful to illustrate the objective of the paper. Fracture length and spacing can play a critical role in an effective hydraulic stimulation treatment, which will have important economic consequences for a potential EGS development in the study area. The authors are aware that the modeling exercise may be uncertain. Further field tests are of upmost importance to decrease the uncertainty and provide more accurate and realistic results.

It is also important to note that carrying out geothermal research in northern communities requires “out of the box” solutions for a first-order assessment of their geothermal potential. Field logistic is extremely challenging in communities not connected by road, where field material (e.g., borehole tools) needs to be shipped by boat or airplane, and the fieldwork timing is reduced due to cost and short field seasons followed by the long winter. If the preliminary data suggests an interesting potential, then it opens the door to further geothermal research.

Another source of bias and uncertainty in the fracture network characterization and modeling exercise comes from the analysts themselves. If other people analyze our datasets and obtain a different interpretation, this may lead to different results, highlighting a reproducibility issue. If different results are obtained, this can become a critical issue as pointed out by Peacock et al. (2019). The data in this paper was analyzed by one single analyst. This by itself is already imposing a subjective bias, and indeed further work needs to be done to combine the interpretation of more analysts and understand the reproducibility of the results. The analyst subjective bias has also been pointed out by Andrews et al. (2019). Nevertheless, this analyst bias was considered not critical to the objective of this paper.

Another important point worth discussing is the numerical simulator used in this study, and the fracture constitutive laws implemented. Joint behavior during hydraulic stimulation is approximated in this study following Goodman/Willis-Richards’ approximation rather than the Barton-Bandis model. This can lead to another source of bias and uncertainty in the modeling results that the authors are aware. In fact, a next step can be the comparison of the results from FRACSIM3D with other numerical simulators to quantify the influence of the numerical simulator on forecasting EGS development and performance.

Despite some bias, uncertainty and limitations related to the difficulty of collecting reliable datasets in remote regions, this research work qualitatively illustrates how the scale effect of fracture geometrical properties can have an important impact on hydroshearing treatment. The results, despite uncertain, suggest that further and more detailed geothermal research focused on the natural fracture network is worthwhile to assess accurately if EGS can be developed in the study area. This research work highlights how neglecting the scale effect of geometrical properties of fractures can have an important economic impact on a deep EGS project.

8 Conclusions

This paper analyzed the natural fractures and structural lineaments in the community of Kuujjuaq (Nunavik, Canada) through different datasets and investigated how the observed variability could impact the network response to hydroshearing through numerical modeling. Despite the scale effect being known, our study illustrates challenges to consider when analyzing fracture data with the goal of generating stochastic DFNs and understand the deep-seated structures.

The results of this work highlight that orientation, size and connectivity are greatly dependent on the scale of observation, and less so for fracture spacing and linear frequency. Outcrop data alone suggests one main fracture set while the combination of satellite image and outcrop data indicates five main families. Furthermore, outcrop data indicates fractures smaller than 30 m while structural lineaments of up to ~ 1,500 m were identified in the satellite image. A power-law distribution provides the best fit for the fracture size. Additionally, core data reveals more widely spaced fractures, and consequently lower linear fracture frequency, than outcrop data. The data also suggests that fracture spacing can be fit to an exponential distribution. Fracture aperture is found to obey a power-law distribution.

The numerical simulations of EGS development and performance suggest that a DFN built exclusively based on satellite image data performs better than a DFN built exclusively based on outcrop data. One reason could be the greater shear displacement experienced by larger fractures, which could be related with their lower shear stiffness. However, a DFN built exclusively based on satellite image data may increase the probability of reservoir thermal short-circuiting. The simulations also illustrate that a DFN built exclusively based on outcrop data may need greater fracture frequency and/or greater stimulation and circulation pressures to achieve similar flow rates to a large fracture size network.

Despite being a case study, the work described here aims at acquiring new knowledge on the geometry of natural fracture networks. This is essential to proceed with exploration steps and evaluate the technical feasibility of EGS in Precambrian basement rocks. This represents a valuable contribution for the geothermal community assessing the potential of geothermal technologies in northern environments and dealing with important data gaps. Deep geothermal exploration boreholes are needed to validate initial estimates based on outcrop analogues. Nevertheless, to raise the interest for geothermal exploration drilling, research needs to be undertaken using data that can be economically collected. This is where this contribution gains its scientific value.

Appendices

Appendix 1

Fracture data was collected in 6 sampling areas around the community of Kuujjuaq (A1 to A6; Fig. 2). In sampling area 1 (A1), the fracture data was collected in 6 outcrops spread around this fracture sampling area and in each outcrop 1 to 3 scanlines were used. The scanlines were placed at an angle to each other to capture more fracture information. This sampling area is characterized by two main fracture directions, E-W and NW–SE, and scattered fractures with orientation N-S and NE-SW (Fig. 

Fig. 11

a Rose diagrams, b fracture length, c linear fracture intensity and d fracture spacing for each outcrop and scanline studied in the sampling area A1. Red lines in the rose diagrams: scanlines, O: outcrop, n: number of fractures

11). The lengths of the fractures range from a minimum of 0.32 m to a maximum of 13 m (Fig. 11). The linear fracture intensity ranges from 0.35 to 4 m⁻¹ (Fig. 11). The average fracture spacing after application of the Terzaghi correction varies from 0.07 to 2.19 m (Fig. 11).

A total of 5 outcrops spread around fracture sampling area 2 (A2) were selected to collect fracture information. One to three scanlines at an angle to each other were used in each outcrop. The main fracture directions are N-S and E-W, and the minor directions are NE-SW and NW–SE (Fig. 

Fig. 12

a Rose diagrams, b fracture length, c linear fracture intensity and d fracture spacing for each outcrop and scanline studied in the sampling area A2. Red lines in the rose diagrams: scanlines, O: outcrop, n: number of fractures

12). The fractures have a length of 0.08 to 5 m (Fig. 12). This sampling area is characterized by a linear fracture intensity varying between 1.11 and 4.89 m⁻¹ (Fig. 12), with fractures spaced by 0.10 to 0.41 m after application of Terzaghi correction (Fig. 12).

Eight outcrops were surveyed in sampling area 3 (A3) by one to four scanlines deployed at an angle to each other. The main fracture direction is E-W, and the minor di-rection is NW–SE, with scattered fractures with orientation N-S and NE-SW (Fig. 

Fig. 13

a Rose diagrams, b fracture length, c linear fracture intensity and d fracture spacing for each outcrop and scanline studied in the sampling area A3. Red lines in the rose diagrams: scanlines, O: outcrop, n: number of fractures

13). The length of the fractures varies between 0.2 and 12 m (Fig. 13), the linear density between 0.48 and 3.54 m⁻¹ (Fig. 13) and the spacing between fractures from 0.001 to 3.46 m after correction (Fig. 13).

The fractures in sampling area 4 (A4) were collected in 3 outcrops and 1 to 4 scan-lines positioned at an angle with each other were used. The main fracture direction identified is E-W and NE-SW, with scattered fractures in the direction N-S and NW–SE (Fig. 

Fig. 14

a Rose diagrams, b fracture length, c linear fracture intensity and d fracture spacing for each outcrop and scanline studied in the sampling area A4. Red lines in the rose diagrams: scanlines, O: outcrop, n: number of fractures

14). The fractures have a minimum length of 0.56 m to a maximum length of 27.78 m (Fig. 14), the linear intensity is analyzed to range between 0.29 and 2.22 m⁻¹ (Fig. 14) and the spacing between 0.05 and 5.46 m, after application of Terzaghi correction (Fig. 14).

Three outcrops were surveyed in sampling area 5 (A5), and one to two scanlines were used to collect fracture information in this sampling area. The scanlines were displaced at an angle to each other. The main fracture direction is N–S and NE–SW with scattered fractures in the direction E-W and NW–SE (Fig. 

Fig. 15

a Rose diagrams, b fracture length, c linear fracture intensity and d fracture spacing for each outcrop and scanline studied in the sampling area A5. Red lines in the rose diagrams: scanlines, O: outcrop, n: number of fractures

15). The length of the fractures in the outcrops analyzed varies between 0.38 and 10 m (Fig. 15). The line-ar fracture intensity is found to range between 1.55 and 1.77 m⁻¹ (Fig. 15), with a corrected fracture spacing varying between 0.07 and 1.60 m (Fig. 15).

Appendix 2

An important step in DFN modeling is the determination of the geometric parameters (i.e., orientation distribution, length distribution, frequency distribution, spacing distribution and aperture distribution) of individual fractures. These parameters serve as the foundation for the simulation.

Orientation distribution

Fisher (Eq. 12; Borradaile 2003), von Mises (Eq. 13; Borradaile 2003; Chabani et al. 2020) and mixture of von Mises (Eq. 14; Chabani et al. 2020) distributions were applied to describe fracture orientation:

$$f\left( a \right) = \frac{k}{4\pi \sinh k}e^{{\left( {k\cos a - \mu } \right)}}$$

(12)

$$f\left( \theta \right) = \frac{1}{{2\pi \left( {I_{0} \left( k \right)} \right)}}e^{{\left( {k\cos \left( {\theta - \mu } \right)} \right)}} ,{\text{ with }} I_{0} \left( k \right) = \frac{1}{2\pi }\mathop \smallint \limits_{0}^{2\pi } e^{k\cos \theta } d\theta {\text{and }}\theta = \frac{a2\pi }{{360}}{ }$$

(13)

$$MvM_{H} \left( \theta \right) = \mathop \sum \limits_{h = 1}^{H} w_{{\text{h}}} f_{{\text{h}}} \left( \theta \right){ }$$

(14)

where k (radians) is a concentration parameter that equals the variance (σ²) of the distribution, μ (radians) is the circular mean (0 ≤ μ < 2π) and a (radians) is the azimuth. I₀ (−) is the modified Bessel function of the first kind and order zero and w (−) is the weight of each component h.

Length/size distribution

The fracture size can follow negative exponential, lognormal, or power-law distributions. The PDF of an exponential distribution is given by (e.g., Zeeb et al. 2013b):

$$f\left( l \right) = 1 - exp\left( { - \left( {\frac{1}{{\mu_{{{\text{length}}}} }}} \right)l} \right)$$

(15)

where l (m) is the fracture length and μ_length (m) is the mean length.

The PDF of the lognormal distribution is (e.g., Zeeb et al. 2013b):

$$f\left( l \right) = \frac{1}{2} + \frac{1}{2}erf\left( {\frac{{\ln \left( l \right) - \mu_{{{\text{length}}}} }}{{\sqrt {2\sigma_{{{\text{length}}}}^{2} } }}} \right)$$

(16)

where erf is the error function and σ_length (m) is the standard deviation.

Of these, the power-law function is widely used to describe the trace length of natural fractures (e.g., Renshaw 1999; Gutierrez and Youn 2015; Lei et al. 2017; Toth et al. 2022):

$$f\left( l \right) = 1 - \left( {\frac{l}{{l_{{{\text{min}}}} }}} \right)^{ - D}$$

(17)

or,

$$N\left( l \right) = c \times l^{ - D}$$

(18)

where N(l) (−) is the number of fractures and c and D are the parameters of the distribution function. The parameter c is equal to:

$$c = \frac{1 - D}{{l_{{{\text{max}}}}^{1 - D} - l_{{{\text{min}}}}^{1 - D} }}$$

(19)

One characteristic of the power-law distribution is that the number of small fractures is significantly larger than the number of longer fractures.

Frequency distribution

Liner fracture intensity (P₁₀; m⁻¹) is defined as the number of fracture intersections per unit length (Zeeb et al. 2013b):

$$P_{10} = \frac{N}{L}$$

(20)

where N (−) is the number of fracture intersections and L (m) is the scanline length. In instances where the first and last fracture do not match the start and end of the scanline, the unmeasured spacing needs to be removed (Sanderson and Peacock 2019):

$$P_{10} = \frac{N - 1}{{d_{N} - d_{1} }}$$

(21)

where d_i (m) is the distance measured from the first fracture d₁ to the last d_N.

Areal fracture density (P₂₀; m⁻²) correspond to the number of fractures per unit area (Zeeb et al. 2013b):

$$P_{20} = \frac{N}{A}$$

(22)

where A (m²) is the area of observation.

Spacing distribution

Terzaghi correction was applied since the fractures sampled were not at right angles with the scanline and that influences the fracture spacing (Zeeb et al. 2013b):

$$S = S_{{{\text{apparent}}}} \times \left| {{\text{cos}}\left( {a - a^{\prime}} \right)} \right|$$

(23)

where a’ (º) is the orientation of the scanline.

Once the spacing is corrected, the mean fracture spacing can be evaluated as (Sanderson and Peacock 2019):

$$S = \frac{L}{N}$$

(24)

where S (m) is the mean spacing. Spacing and intensity are related as (Zeeb et al. 2013b):

$$S = \frac{1}{{P_{10} }}$$

(25)

The spacing parameter is often described by the coefficient of variation (e.g., Gillespie et al. 1993; Sanderson and Peacock 2019):

$$CV = \frac{{\sigma_{{{\text{spacing}}}} }}{{\mu_{{{\text{spacing}}}} }}$$

(26)

where CV (−) is the coefficient of variation and σ_spacing (m) and μ_spacing (m) are the standard deviation and mean spacing, respectively. The coefficient of variation provides an indication of the fracture distribution and clustering (Table

Table 4 Distribution and clustering of fractures according to the coefficient of variation (e.g., Gillespie et al., 1993; Manzocchi 2002)

4).

Aperture distribution

The aperture (w; m) of the fractures evaluated on the outcrops can follow lognormal or power-law distributions. Both were applied to the data to evaluate which has the best fit:

$$f\left( w \right) = \frac{1}{2} + \frac{1}{2}erf\left( {\frac{{\ln \left( w \right) - \mu_{{{\text{aperture}}}} }}{{\sqrt {2\sigma_{{{\text{aperture}}}}^{2} } }}} \right)$$

(27)

$$f\left( w \right) = 1 - \left( {\frac{w}{{w_{{{\text{min}}}} }}} \right)^{ - D}$$

(28)

Appendix 3

For each geometric parameter, several probability distribution functions (PDFs) were applied to assess which provides the best-fit to the data sampled. Statistical tests such as the root mean squared error (RMSE; Eq. 29) and the sum of squared errors (SSE; Eq. 30) were used to evaluate the goodness of the fit (e.g., Zeeb et al. 2013b):

$$RMSE = \sqrt {\frac{{\mathop \sum \nolimits_{i = 1}^{N} \left( {P_{i} - O_{1} } \right)^{2} }}{N}}$$

(29)

$$SSE = \mathop \sum \limits_{i = 1}^{N} \left( {P_{i} - O_{1} } \right)^{2}$$

(30)

where O is the sampled/observed value and P is the estimated/predicted value.