Coupled Evaluation of DFN Models and Well Hydraulics for Improved Estimation of Hydraulic Fracture Apertures

1. Introduction

A deeper understanding of the hydrodynamic behaviour of fractured rock masses is essential in many applied geology fields (e.g., [1,2,3]). In these formations, fluids can move and accumulate within the complex network of fractures rather than in pores. In purely fractured rock masses, matrix porosity is minimal; however, in other situations, both intergranular and fractured porosity may coexist in varying proportions. While fractured sandstones often exhibit a dual-porosity system, most igneous and metamorphic rocks are generally considered fractured. The gradual transition between fully porous and fully fractured states is described by Nelson [4]. The network of discontinuities typically forms through structural deformation; however, atectonic processes, such as cooling in volcanic rocks or the development of joint systems within any rock body (e.g., [5], can also induce fracture formation. Regardless of their origin, many fracture features vital to their hydrodynamic role may undergo significant modifications after initially forming. Key processes include water-rock interactions, dissolution of the host rock, precipitation of new mineral phases filling veins [6], and hydration of primary minerals within the host rock [7]. Moreover, the geometry of individual fractures, and consequently the entire fracture network, can change with depth (either subsidence or uplift; [8]), affected by the anisotropy of the recent stress field [9,10,11], or even by fluid overpressure within the reservoir [12]. All these processes simultaneously influence the geometry of both the fractures and their network. Laboratory experiments [13] on limestone specimens demonstrate that percolating water can cause the dissolution of the host rock in certain sections of a fracture system, while the precipitation of new vein-filling minerals in other sections results in ongoing geometric changes. In the case of an anisotropic recent stress field, the tendency for fractures to close highly depends on their orientation relative to the principal stress axes. Fractures aligned with σ1 tend to open, whereas those perpendicular tend to close (e.g., [14].

Due to the complex interaction between tectonic processes and fluid-rock interactions, fracture surfaces are usually rough. This roughness can significantly affect fluid seepage behaviour within jointed rock masses; the rough surface of the fracture wall may reduce the aperture through which fluid actively moves [15,16]. In fact, if a fracture wall surface is uneven, the hypothesis of the parallel plate model for describing its hydrodynamic behaviour becomes invalid [17]. Several practical methods exist to address this issue; for calculating reliable fracture apertures, correction factors such as the Joint Roughness Coefficient, Fractal Dimension, Hurst number, and Asperity Height are commonly employed [18,19,20,21]. Changes in the volume of individual fractures influence the fundamental hydrodynamic parameters, porosity and permeability, which also depend on the combined effects of primary structural evolution and post-tectonic modifications. Since all these changes have minimal impact on fracture length [22], variations in spatial extent are best addressed by modifying fracture aperture.

However, determining not only the measurement but also the correct definition of a fracture’s aperture is challenging ([10,23], and references therein). Since all fractures are somewhat curved and have rough surfaces [24], a simple definition could be the average distance between the opposite rock walls (physical aperture). This measurement can be taken on the surface at any scale (outcrop, hand specimen, microscopy) with sufficient accuracy [25]. However, the measured values are irrelevant under subsurface reservoir conditions because of fracture closure caused by lithostatic pressure increasing with depth and the influence of recent stresses. Some well-log methods (acoustic, BHTV, etc.) may provide a reliable approximation of the physical aperture underground, above a specific fracture size [26,27].

Instead of using a physical aperture, paleo-apertures can be employed for modelling, defined as the width of mineralised veins within rocks [28]. Measuring this parameter is straightforward on the surface, as the extent of the veins may not change with surface uplift and has remained constant since their formation. Therefore, they retain the original apertures. However, they should not be used to represent recent reservoir conditions. Another common approach for hydrodynamic modelling is the constant aperture method [10,29]. In this method, each fracture is assigned the same width regardless of its lateral size.

From a hydrodynamic perspective, the most critical parameter is hydraulic aperture, which, by definition, is the aperture of an ideal, smooth-walled conduit that would produce the same flow rate under a given pressure gradient as the actual fracture [27]. Although this value appears to function as a geometric parameter, it cannot be measured directly through geometric means. To model the behaviour of a fractured reservoir more accurately, valid hydraulic apertures should be estimated based on current reservoir conditions rather than any other definition of fracture aperture [30]. Since hydraulic aperture is not solely a geometric variable, well-hydraulic data should be included in the modelling process to estimate it. These essential raw data typically come from pumping tests, which provide transmissivity (hydraulic conductivity) values for well-defined sections of the reservoir.

A standard method for modelling fractured rock bodies is to represent the fracture network as a geometric pattern. This approach utilises a well-defined set of geometric parameters [31], including fracture length distribution, the spatial location of fracture centres, and fracture orientation (dip, dip angle), among others (spacing, intersection, etc.). These data provide a solid foundation for most fracture network modelling algorithms. The most widely used methods are DFN (discrete fracture network) techniques [32,33,34], which are based on different mathematical frameworks. In this project, the RepSim code [22,31,35,36,37] was employed.

This approach assumes that length data follow a power-law distribution [38] and that the spatial positions of the fracture midpoints can be approximated as a fractal-like pattern. The mathematical basis of the software code is explained in detail by M. Tóth et al. [39]. Using a DFN simulation algorithm, the geometry of the fracture network can be reconstructed in three dimensions. Since many of these algorithms are stochastic, multiple models with the same probability can be simulated and assessed. Such models enable answers to key questions about fractured rock bodies, such as the connectivity between individual fractures and the sizes and spatial arrangements of communicating sub-systems.

Calculating reliable porosity and permeability data depends not only on a realistic 3D fracture network model but also on selecting appropriate aperture values for individual fractures. Most observations indicate a steady increase in aperture with fracture length; therefore, when length data are available, the aperture can be determined directly. Some methods employ a linear model, while others favour an exponential function for this approximation. Although this choice significantly influences the initial aperture values, it cannot account for the effects of complex post-tectonic processes that can markedly alter the aperture. As late-stage processes alter fracture apertures, the simple relationship between length and aperture becomes unreliable, making the accurate estimation of aperture data more challenging. Although such data are essential for accurately assessing the porosity and permeability of a fractured rock body, directly measuring or even estimating typical hydraulic aperture values remains a challenge.

Conversely, transmissivity can be determined through a standard pumping test conducted at suitable intervals in any well that penetrates a fractured rock formation. These data can be converted into intrinsic permeability values, which can also be estimated based on any simulated 3D fracture network [40,41]. Since the measured and estimated permeability data are independent, and the latter depends solely on the chosen aperture values, comparing the measured and estimated permeabilities can help calibrate apertures under reservoir conditions.

This paper aims to demonstrate how to accurately calculate hydraulic aperture data for fractured rock bodies by integrating measured transmissivities into DFN models, and to utilise these data more effectively to better understand the hydrodynamic behaviour of a fractured reservoir.

2. Methods

2.1. DFN Modelling

Most fracture network simulation software employs a DFN algorithm and utilises quantifiable geometric parameters to model fractures. Among these, trace length distribution, fracture orientation data, and a parameter that evaluates the spatial density of fractures are the most vital. Length data are generally assumed to follow a power-law distribution so that N(L) = F*L^-E, where N(L) is the number of fractures of length L, while E and F are constants [38,42]. This relationship enables the estimation of the number of fractures of any length, including those outside the measured length range.

The spatial position of a single fracture can be most accurately estimated using the fractal dimension of its midpoint pattern in 3D [43,44]. This parameter can be assessed using the box-counting method, in which the fracture midpoint pattern is first overlaid with a regular grid of square boxes. Then, the number of boxes (N(r)) needed to cover the pattern is counted. The fractal dimension (D) is determined by how this number increases as the grid becomes finer: N(r) ~ r^-D. The key geometric parameters, D and E, can be obtained by analysing any 2-dimensional section of a fractured rock body or by using 1D sections well-logs [36,37,45].

Using these parameters, along with a set of dip (α) and dip angle (β) values, the main steps of modelling in the RepSim package are as follows:

1)

Definition of the study area and a regular grid network.

2)

Measurement of the primary geometric parameters (D, E, F, α, β) for several localities (grid cells) using either outcrops or well-logs.

3)

Using any interpolation algorithm (e.g., nearest neighbour, inverse distance, ordinary kriging, etc.) to calculate reliable values for all the above parameters for all previously unknown grid cells.

4)

Based on the primary parameters, fracture networks can be stochastically generated for the entire study area. As a result, any number of realisations (fracture networks in 3D) with identical probability can be created. All these networks mirror the geometry of the original fracture pattern in the study area. For details on the discrete fracture network generation process using the RepSim code, see [22,31,35,45].

Finally, the aperture can be calculated for any single fracture using a common deterministic relationship. Typically, a = A*L is used [43,46,47,48], where a is the aperture and A is a constant. From now on, A is referred to as the aperture coefficient. Although the linear model is widely adopted, Olson [49] notes that the relationship between length and aperture is often nonlinear for fractures composed of multiple segments and for strata-bound fractures. An aperture calibration algorithm similar to the one discussed in the following chapters could also be developed for this approximation.

2.2. Aperture Calibration

Having adequate geometric data on the fracture network offers a strong foundation for simulating a DFN model for any part of the rock mass. By using such a model and choosing a reliable aperture coefficient, both fractured porosity and the intrinsic permeability tensor can be determined for any cell within the modelled rock mass. The algorithms and formulae for estimating these parameters follow those of [35].

Fractured porosity is defined as

$$\mathrm{\Phi }={\displaystyle \frac{Vf}{V}}.$$

(1)

In the case of cubic cells, $V=r^3$, while the total volume of the fractures within a particular cube (V_f) can be accurately approximated by the lower Riemann sum, that is,

$$Vf={\displaystyle \sum _{i=1}^n{\displaystyle \sum _j\frac{l_{ij}\cdot a_{ij}\cdot r}{n},}}$$

(2)

where j is the number of fractures intersected by the ith slice in the Riemann algorithm (Figure 1), and porosity appears as

$$\mathrm{\Phi }={\displaystyle \frac{1}{\mathrm{n}*{\mathrm{r}}^2}}{\displaystyle \sum _{\mathrm{i}=1}^n{\displaystyle \sum _{\mathrm{j}}{\mathrm{l}}_{\mathrm{i}\mathrm{j}}\cdot {\mathrm{a}}_{\mathrm{i}\mathrm{j}}}}$$

(3)

The permeability of a fractured rock mass can be described by a 3×3 permeability tensor. In the RepSim code, it is computed using Oda’s [50] algorithm. Therefore, under Darcy’s law,

$${\mathrm{v}}_{\mathrm{i}}={\displaystyle \frac{\mathrm{g}}{\mathrm{\mu }}}\cdot \mathrm{\rho }\cdot {\mathrm{k}}_{\mathrm{i},\mathrm{j}}\cdot {\mathrm{J}}_{\mathrm{i}}$$

(4)

where v is the specific flow rate, μ is the dynamic viscosity, ρ is the fluid density, and J is the hydraulic gradient. Conversely, since fluid can only percolate along fractures within a given volume,

$${\mathrm{v}}_{\mathrm{i}}={\displaystyle \frac{1}{\mathrm{V}}}\cdot {\displaystyle \int {{\mathrm{v}}_{\mathrm{i}}}^{\mathrm{f}}\cdot \mathrm{d}{\mathrm{V}}^{\mathrm{f}}}$$

(5)

where ${\mathrm{v}}^{\mathrm{f}}$ is the flow velocity in a discrete fracture. This formula is approximated ad libitum by

$${\mathrm{v}}_{\mathrm{i}}={\displaystyle \frac{1}{\mathrm{V}}}\cdot {\displaystyle \sum _{\mathrm{f}}{{\mathrm{v}}_{\mathrm{i}}}^{\mathrm{f}}\cdot {\mathrm{V}}^{\mathrm{f}}}$$

(6)

Under the cubic law, where assuming laminar flow within a fracture (parallel plate model, [51,52], the specific flow rate is proportional to the square of the fracture aperture, and

$${{\mathrm{v}}_{\mathrm{i}}}^{\mathrm{f}}={\displaystyle \frac{1}{12}}\cdot {\displaystyle \frac{\mathrm{g}}{\mathrm{\mu }}}\cdot \mathrm{\rho }\cdot {\mathrm{a}}^2\cdot {{\mathrm{J}}_{\mathrm{i}}}^{\mathrm{f}},$$

(7)

where ${{\mathrm{J}}_{\mathrm{i}}}^{\mathrm{f}}$ is the i-th component of J projected onto the f fracture, that is, as

$$J^f=J-(n\cdot J)\cdot n\mathrm{and}$$

(8)

$${J_i}^f={\displaystyle \sum _j({\delta }_{ij}-n_i{n}_j)\cdot J_j},$$

(9)

where ${\mathrm{\delta }}_{\mathrm{i}\mathrm{j}}$ is the Kronecker delta symbol. Thus, finally, comparing (4) and (6) according to Oda [50],

$${\mathrm{k}}_{\mathrm{i},\mathrm{j}}={\displaystyle \frac{1}{12}}\cdot ({\mathrm{P}}_{\mathrm{k}\mathrm{k}}\cdot {\mathrm{\delta }}_{\mathrm{i}\mathrm{j}}-{\mathrm{P}}_{\mathrm{i}\mathrm{j}})$$

(10)

and under the discretisation solution of Koike & Ichikawa [53], considering that in the case of cubic cells, V = r³,

$$P_{ij}={\displaystyle \frac{1}{r^3}}\cdot {\displaystyle \sum _f{a}^3\cdot l\cdot n_i\cdot n_j}$$

(11)

Finally, using the lower Riemann sum approximation

$$P_{ij}={\displaystyle \frac{1}{k\cdot r^3}}\cdot {\displaystyle \sum _k{\displaystyle \sum _f{a}^3\cdot l\cdot n_i\cdot n_j}}$$

(12)

and

$$P_{kk}=P_{11}+P_{22}+P_{33},$$

(13)

where n_i and n_j are the normal vector projections of the given fracture onto the particular axes.

It follows from these formulae that in any DFN model (with fixed D, E, F, α, and β parameters), the values of both the calculated porosity and permeability depend solely on the aperture, and ultimately on the A aperture coefficient.

Hydraulic data indicating the hydrodynamic behaviour of a fractured system can be obtained through various laboratory and field methods. By analysing the raw data from in situ well-hydraulic measurements, such as different pumping tests, hydrodynamic parameters (including hydraulic conductivity and transmissivity) are estimated to characterise specific sections of the reservoir. Relevant details on measurement procedures and evaluation methods can be found in the literature [54,55], though these are outside the scope of the present paper and are not discussed here. Hydraulic conductivity (K) and transmissivity (T) are directly proportional, with T = K × d, where d is the thickness of the saturated interval. Additionally, hydraulic conductivity can be converted into permeability (k) using

$$k={\displaystyle \frac{K·\mu }{\rho ·g}},$$

(14)

where μ is the dynamic viscosity of the fluid, ρ is the density, while g stands for gravitational acceleration.

A permeability value derived from measured hydrodynamic parameters relates to a specific segment of the rock mass. If the geometric parameters of the fracture network (D, E, F, α, and β) are available, it is possible to generate a DFN model valid for that same segment. Using the fracture aperture coefficient (A), the intrinsic permeability tensor for this segment can be calculated (formulae 4–14), enabling the observed and estimated permeabilities to be compared. While the observed permeability depends on the well-hydraulic measurement and its evaluation method, the estimated permeability relies solely on the A parameter. Therefore, comparing observed and estimated permeability values can help calibrate the aperture coefficient (Figure 2) to align them as closely as possible. During calibration, the A coefficient should be identified for which log(k_obs)/log(k_est) = 1.

3. Case Study

3.1. Geological Background

To introduce the approach described above, the Mórágy Granite Formation (MGF) in southwest Hungary was selected (Figure 3). This formation is well known for its geological characteristics and was chosen as the repository for low- and medium-level radioactive waste in Hungary. During site preparation for the repository, over 20 boreholes were drilled through the granitoid massif, achieving nearly 100% core recovery. Additionally, two underground access tunnels were excavated. The results of the petrological, structural geological, and hydrogeological evaluations of the rock body are summarised in several previous papers [56,57,58,59] and are briefly presented hereafter.

The intrusive mass consists of two types of Lower Carboniferous granitoid rocks (Figure 3); the porphyritic monzogranite is interbedded with a more mafic rock of monzonitic composition. Recent models suggest igneous evolution through a magma-mixing process. Two successive ductile events have influenced the early structural history of the rock body. In addition to the overall NE-SW-striking magmatic foliation, steeply foliated mylonitic zones with the same strike are observed [56]. The brittle deformation history of the massif is complex, as evidenced by multiple deformation phases in both outcrop and borehole core samples [60]. The main large-scale deformation zones exhibit two typical orientations: the dominant NE-SW set and another perpendicular set (NW-SE). During 3D fracture network modelling, the geometric parameters (D, E, F, α, β) of the fracture system were measured and evaluated at every 25-metre interval across 14 wells in the study area [35]. Consequently, near-well DFN models were generated for each well. As an example, calculated D values and the final fracture network model for well Üh-2 are shown in Figure 4. These models reveal a highly heterogeneous fracture system within the granitoid body, with more- and less-fractured segments along each well. Regarding hydrogeological implications, two hydraulic domains can be distinguished within the fractured granitoid body. The less transmissive blocks (LTB) are separated by more transmissive zones (MTZ), which follow the major tectonic lines in NE-SW and NW-SE directions [59].

Studying the hydrodynamic behaviour of single wells, Benedek & Dankó [58] highlighted the coexistence of small-scale hydraulic head scattering and large hydraulic head jumps along individual boreholes. The previous DFN models (Figure 6 in [35]) indicate that the hydraulic head tends to change abruptly at specific depth horizons, where a connected fracture network could have formed. These model results seem to confirm the conclusion of Balla et al. [61], who suggest that abrupt head jumps are caused by highly altered fault core zones rather than a sparse fracture network. Consequently, mineralisation from water-rock interactions could have significantly altered the original apertures of individual fractures.

A typical aperture coefficient for fractured granitoid rock was determined from hand specimens under the microscope, yielding A = a/L ~ 2.7×10⁻² on average [35]. Using the above DFN models with these measured fracture aperture data, the eigenvalues of the intrinsic permeability tensor average 2.34×10⁻¹⁴, 1.89×10⁻¹⁴, and 1.22×10⁻¹⁴ m². These values are three orders of magnitude larger in the most heavily fractured zones, at 1.71×10⁻¹¹, 1.62×10⁻¹¹, and 1.28×10⁻¹¹ m². Employing the same models and aperture data, the average porosity is 1.62%, with a maximum of approximately 6.00% [35]. However, all these results are based on raw physical aperture measurements taken on the surface and may differ significantly from the actual values. This calculation does not account for corrections due to depth, stress field anisotropy, or the influence of post-tectonic alteration processes.

3.2. Results of the Aperture Calibration

Aperture Coefficient Values

As part of a detailed investigation into the radioactive waste repository, in situ hydrodynamic measurements were performed in numerous wells. Data were collected from a total of 238 depth segments across ten wells and analysed (Figure 3). Transmissivities were measured over approximately 10 m long intervals and used to calculate permeabilities. The measured permeabilities vary by 3–4 orders of magnitude across all studied wells, ranging from a minimum of 10-16 to 10-18 m² to a maximum of 10^-14 to 10^-12 m², supporting the heterogeneous structure of the rock mass (Table 1). Harmonic averages, which represent typical values, are around 10^-15 to 10^-17 m², that is, 1–3 orders of magnitude lower than those calculated using physical aperture data [35].

All fracture network parameters (D, E, F, α, and β) were determined using BHTV and core scanner data [62,63] and are available for all sections used for hydraulic measurement in each studied well [35]. Based on this, DFN models were independently simulated for all 238 intervals, employing the algorithm described above. A permeability tensor was calculated from these fracture network models, assuming A=0.1, 0.08, 0.06, 0.04, and 0.02, respectively, for each studied interval. The diagonal elements of the intrinsic permeability tensor vary in all cases, confirming the anisotropic geometry of the fracture network. However, since field measurements do not provide information about the orientation of permeability (transmissivity), the eigenvalues of the tensor were averaged to enable comparison between measured and observed data. By fitting a logarithmic curve to the points on the observed (k_obs) and estimated (k_est) permeabilities, the aperture coefficient (A) for each interval can be computed, where log(k_obs)/log(k_est) = 1 (Figure 5). The A-values calculated fluctuate within a broad range across all wells. For example, in a single well (Üh-2), the A-values range from 1.4E-4 to 3.3E-3, indicating that a fracture with a diameter of 1 m could have a hydraulic aperture (a = A*L) ranging from 0.14 mm to 3.3 mm.

4. Discussion

As the above computations show, aperture data vary across several orders of magnitude in all wells, likely reflecting the different intensities of post-tectonic aperture-altering processes such as cementation, dissolution, and closure. Without a detailed petrographic analysis, it is difficult to assess the relative significance of these independent effects. However, the spatial pattern along a single well (Üh-2) suggests that the combined influence of these processes defines intervals with typical behaviours (Figure 6a). In sections where A tends to increase, dissolution may be a key process, or late cementation effects are less important. Conversely, zones with unusually low A-values probably represent cemented intervals within the rock formation.

4.1. Calculated Porosities

Such differences should also lead to variations in hydrodynamic characteristics at different parts of the well. To understand the hydrodynamic behaviour of the fracture system, the fractured porosity should be calculated using the fracture model and the derived A parameter for each interval (see formulae 1–3). Based on this data, a porosity log can be created along the wells using the calibrated aperture coefficients (see Figure 6b). This porosity log differs from those estimated using an average A-value or obtained by evaluating well-log geophysical signals. It more accurately reflects the actual hydrodynamic behaviour of the rock mass and aligns more closely with measured transmissivity (permeability) data. Finally, measured permeability and calculated porosity data can be plotted on a traditional porosity-permeability cross-plot and analysed in any conventional manner (e.g., [64] and references therein; [65,66].

Figure 6. a) Spatial pattern of the aperture coefficient (A) along the Üh-2 well indicates that the outcome of post-tectonic processes determines intervals of typical behaviours. b) A porosity log can be sketched along the well using the calibrated aperture coefficients. c) The rock body penetrated can be divided into four major zones hydrodynamically based on the FZI values.

Figure 6. a) Spatial pattern of the aperture coefficient (A) along the Üh-2 well indicates that the outcome of post-tectonic processes determines intervals of typical behaviours. b) A porosity log can be sketched along the well using the calibrated aperture coefficients. c) The rock body penetrated can be divided into four major zones hydrodynamically based on the FZI values.

4.2. Evaluation of the Approach Introduced

Numerous methods exist to estimate how aperture varies during different WRI processes. Most experiments examine the behaviour of individual fractures. Dissolution and precipitation studies in carbonate rock samples demonstrate that, while dissolution may dominate in one part of a fracture, mineral precipitation can occur elsewhere. These small-scale variations can be measured using parameters such as the roughness coefficient [67], which are obtained through precise laboratory tests. Naturally, these results reflect the combined influence of several concurrent processes (dissolution, vein cementation) that can alter the shape of the vein wall, even within a single fracture. There are also formulae to estimate how individual fractures close with increasing depth due to lithostatic pressure [8]. Other studies aim to approximate the influence of recent stress-field anisotropy on fracture closure and/or opening [11]. However, as the angle between the principal stress axes and a fracture changes, the overall effect of stress anisotropy cannot be reliably assessed at the reservoir scale. Therefore, all these estimates of aperture change are valid on a small scale but become highly uncertain at the reservoir level, where WRI processes and pressure effects interact. Additionally, their influence can vary significantly across different points within a real rock mass. Consequently, the current approach does not attempt to understand each factor independently at a small scale. Instead, it focuses on characterising their combined effect mathematically.

In the current approach, both porosity and permeability are connected to the volume of the rock body, as measured through well hydraulics. Assuming this size is sufficiently close to the representative elementary volume (REV, [68], these hydrodynamic parameters accurately characterise the fracture system. Therefore, neither of them relates to small samples analysed under laboratory conditions.

4.3. The Case Study

On the poro-perm plot of the Üh-2 well (Figure 7), two distinct trends can be identified, both characterised by >0.9 correlation coefficients using the k = α*ϕ^β regression function, where k denotes permeability, ϕ is effective porosity, and α and β are constants.

The relationship between porosity and permeability is crucial for understanding fluid flow through porous media. It is well established that higher porosity generally correlates with higher permeability. This relationship is often expressed mathematically using empirical equations, typically in the form of k = α*ϕ^β, where α and β are parameters determined experimentally [69,70]. Power-law relations have also been observed across various lithologies, including sandstones and limestones [71,72]. However, it is crucial to recognise that the specific values of α and β can vary considerably depending on material properties and stress conditions [73,74]. Recent research indicates that although porosity is key to predicting permeability, effective porosity (the proportion of pore space that actively enables fluid flow) generally provides a more accurate prediction than total porosity, particularly in complex rock formations [75,76]. Furthermore, the stress field significantly influences both porosity and permeability. Under different stress conditions, pore structures within materials can change in ways that either enhance or reduce permeability [77,78]. For example, the concept of effective stress is crucial for understanding how external pressures influence the permeability-porosity relationship, particularly in shale gas reservoirs and similar geological settings [78,79]. As stress increases, it often results in pore compression, leading to decreases in both porosity and permeability [77,79].

The complexity of the relationship between porosity and permeability is also evident across various rock types. For instance, volcanic rocks show diverse permeability behaviours at different porosity levels, indicating that microstructural heterogeneity can significantly influence fluid flow paths [80,81]. Similar complexity is observed in various hard rock bodies, where fractures can significantly impact flow pathways [72,82].

In conclusion, although a direct correlation usually exists between porosity and permeability, this relationship is influenced by multiple factors, including rock type, stress conditions, and the concept of effective porosity. The importance of this relationship extends beyond theoretical models, as it holds practical significance in reservoir modelling based on measured field data.

In addition to the two trends identified in Figure 7, significantly different permeabilities are associated with the same porosity values, possibly due to varying intensities of post-deformational processes. Furthermore, segments characterised by the two distinct trends form continuous intervals along the well, suggesting that the intensities of these post-tectonic processes also vary spatially; there must be zones that are highly cemented and others mainly characterised by open fractures.

Using effective porosity and permeability data together, their variation along the well can be characterised by several parameters. Among these, the FZI (Flow Zone Indicator) is one of the most indicative of reservoir quality. Amaefule et al. [83] introduced a method to calculate the Flow Zone Indicator so that FZI = sqrt(k/ϕ)/ϕ. The lower this parameter is, the more closed the fractures are. Considering the Üh-2 well discussed above, the FZI data and trends in the poro-perm plot (Figure 7) suggest that the rock body penetrated can be divided into two main types hydrodynamically (Figure 6c). The low FZI values and the coinciding low β exponent of the fitted power function trend of the upper approximately 50 metres imply cemented fractures. This interval is followed by a 20 m short zone with high FZI values and a large β value. The coincident good reservoir quality is likely due to the presence of a hydrodynamically active fault zone. The lowermost 200 metres of the well show a continuous increase in FZI values from approximately 600 up to about 3000 μm. In this segment, the flow characteristics of the fracture system improve with depth. Whether this unusual phenomenon results from hindered vein cementation or from fractures remaining open due to the recent stress field cannot be determined without a detailed geophysical and petrological study. Nevertheless, this 200 m long interval can be divided into two sections based on the poro-perm diagram (Figure 7), indicating that a more permeable trend characterises the lower, approximately 100 m section.

The data from another nine analysed wells can be interpreted similarly. The A-parameter (aperture coefficient) generally varies by 1–2 orders of magnitude within each well, indicating significant differences in the intensities of post-tectonic processes along the wells (Figure 8). High A values clearly indicate segments with more open fractures, while low A values suggest more closed fractures.

Porosity logs further support this, showing sections with higher and lower porosity along each well (Figure 9). On poro-perm diagrams, these sections usually outline independent trends of the relation (k = α*ϕ^β) between porosity and permeability (e.g., Figure 10 and Figure 11). Such trends highlight segments with similar hydrodynamic behaviours. It is important to note that in each case, data points exhibiting similar hydrodynamic behaviour are not randomly distributed across the wells; instead, they appear side by side and define specific depth intervals. In these intervals, the effectiveness of petrological and/or structural processes that influenced the aperture must have been characteristic. These sections can be described as hydraulic flow units (HFU; [83,84], which are best identified along the wells using their FZI logs [85] (Figure 12).

The relationship between the Flow Zone Indicator (FZI) values and Hydraulic Flow Units (HFUs) in fractured rock masses is a key concept in reservoir geology. The FZI characterises rock properties related to fluid flow, which is essential for understanding and predicting reservoir performance. HFUs are defined as volumes of rock where the geological and petrophysical properties that influence fluid movement are internally consistent [86]. This uniformity is crucial because it determines how fluid interacts with the rock matrix. As the FZI is a composite index that combines the rock’s porosity and permeability, it offers insight into the flow characteristics within a given HFU. The FZI value reflects the hydraulic properties of the rock; for example, a higher FZI value generally indicates better reservoir quality and higher permeability, thereby improving fluid flow potential.

Research has shown that different rock types and geological conditions influence FZI values and, consequently, the classification of HFUs within fractured rock masses [36,86]. Several studies demonstrate the applicability of FZI in predicting HFUs. For example, [87] found that incorporating FZI calculations significantly improved the identification of HFUs in the Bredasdorp Basin, indicating that consistent FZI values reflect similar hydraulic flow characteristics [87,88]. Likewise, [89] emphasised the usefulness of FZI in enhancing the accuracy of permeability predictions based on rock typing in carbonate formations.

Considering all calculated data, the FZI values of the Mórágy granite are very high compared with those typical of intergranular-porosity reservoirs (mainly sandstones). However, they fall within the common range for other fractured rock masses, such as carbonates [90,91]. Nonetheless, the specific FZI values are inconsistent, varying by two orders of magnitude across the studied wells, with an FZI of approximately 40 μm for Üh-4. In contrast, the lowermost part of the Üh-22 well reaches a maximum of 4000 μm.

Based on the typical vertical pattern with well-defined sections of markedly different FZI values, all studied wells can be divided into numerous HFUs. In the case of Üh-3 (Figure 12a), FZI~500 m is consistent throughout the well, with two zones of unusually low FZI values. The sudden appearance of these zones indicates the presence of fault zones hosting highly mineralised veins. All porosity and permeability data follow the same pattern for Üh-4 (Figure 10), suggesting a similar post-tectonic evolution of the fracture system around it. Low FZI values are consistent throughout the crosscut segment, indicating uniform hydraulic behaviour of the surrounding rock mass. In the case of the Üh-5 well, three distinct trends can be identified in the poro-perm plot (Figure 11). Samples belonging to each of these groups represent well-defined depth intervals along the well, with FZ indicator values decreasing steadily downward (Figure 12c). Such behaviour reflects fracture closure with increasing depth, due to rising lithostatic pressure or enhanced significance for vein mineralisation.

In the fracture system of the Üh-22 well, the measured permeabilities vary across a range of four orders of magnitude (Table 1, Figure 13a). Simultaneously, the estimated porosities range between 0.05% and 1.2% (Figure 13b). The coinciding poro-perm plot (Figure 14) indicates three distinct trends, all of which correspond to separate depth intervals (HFUs). A continuously decreasing FZI trend is characteristic of the well (Figure 13 c), reaching its minimum at approximately -50 m (a.s.l.), where very low aperture coefficients (~5E-4) are typical, likely due to intensive vein mineralisation. Further down, the sharply increasing FZI values (from ~300 μm to ~3000 μm) denote a new hydraulic flow unit at the lowest studied interval. This pattern is not evident on either the permeability or porosity logs (Figure 13a, b), but it reliably correlates with the observed hydrodynamic behaviour of the well. Comparing the FZI log with the hydraulic head profile of the same well [58] Figure 13d) reveals an apparent coincidence: near-uniform head values extend down to about -140 m (a.s.l.), followed by an abrupt drop in head at the bottom 70 m of the well, where anomalously open fractures in the new HFU can be inferred. This pattern supports the model proposed by [61], which suggests that the lowermost section of the well is isolated from the upper part by an intensively altered and mutually fractured horizon, rather than by a sparse fracture network. Benedek et al. [92] also emphasise that the compartmentalised behaviour of the granitoid reservoir often results from extensive secondary mineralisation of specific fracture zones.

The Üh-25 (Figure 12d) well exhibits behaviour similar to that of Üh-3. Its profile, characterised by a typical FZI value, is interrupted by a roughly 20 m wide zone of low values. This zone can likely be interpreted as a highly mineralised fault zone. The behaviour of Üh-26 (Figure 12e) resembles that of Üh-2 and Üh-22, with a nearly constant FZI trend that decreases downward and two zones of higher values, indicating high-permeability fault zones. An apparent bimodality is evident in the case of Üh-27 (Figure 12f). High FZI values dominate the upper part of the well, whereas lower values are typical as the well extends downwards. Given that the well penetrates the granitoid massif at the boundary between two lithologically distinct realms (Figure 3), the vertical pattern appears to reflect these lithological and rheological changes.

A previous evaluation [35] showed that the two lithologies could be characterised by significantly different fracture network parameters (D, E) and connectivity relationships. Consequently, the hydrodynamic behaviour of the fracture network in the monzogranite and monzonite-dominated blocks should also differ. Along with the Üh-28 well (Figure 12g), a continuous decrease in FZI values is typical, with no sudden change, suggesting that the well probably did not intersect any major fault zone. However, the data define two HFUs, as the values reach a threshold at approximately 100 m (a.s.l.), and the data follow a different trend in the poro-perm plot (Figure 15). The FZI pattern exhibited by Üh-29 (Figure 12h) is similar to those of Üh-2, Üh-22, and Üh-26. A consistent decreasing FZI trend is observed throughout the entire well, with two roughly 20 m wide intervals where the values jump abruptly. The very high FZI values in these intervals suggest hydraulically active fault zones at these depths.

Considering all studied wells, the fractured Mórágy granite can be characterised by the simultaneous presence of numerous hydraulic flow units, with their boundaries identifiable using FZI logs in each case. Most wells penetrating the granitoid body exhibit homogeneous or slightly decreasing FZI values downwards, interrupted by sections sometimes 10 m wide, where an abrupt change in FZI values is typical. These anomalous zones also define independent trends on the poro-perm plot. While the continuous decline in FZI may reflect fracture closure due to increasing lithostatic pressure, the jumps along the logs indicate discontinuous changes. In the case of Üh-27, such variation in FZI values can be interpreted as resulting from two alternating granitoid lithologies. Nevertheless, all other cases likely indicate the presence of a fault zone. In agreement with previous structural geological and hydrogeological studies [58], some fault zones in the area are highly mineralised and behave as impermeable barriers, while others serve as efficient groundwater conduits. The vertical profiles shown by Üh-2, -22, -26, and -29 align well with the widely accepted hydrodynamic model of the region, in which less-transmissive blocks (fractured rock bodies) coexist with more-transmissive zones (fault zones). Conversely, the vertical patterns of Üh-3 and -25 exemplify the presence of intensely mineralised fault zones.

A more detailed analysis of the aperture, effective porosity, and FZI logs of all studied wells is beyond the scope of this paper, as it would require a more specific geological background. The limited examples provided, however, demonstrate the applicability of the introduced algorithm for estimating the aperture coefficient, fractured porosity, and the Flow Zone Indicator by simultaneously using measured transmissivity data and DFN models.

5. Conclusions

Overall, the approach described above offers many advantages. Calculating reliable hydraulic aperture coefficients from fracture geometry data and measured transmissivities provides crucial insights into how fractures open and close in the studied wells. The fracture aperture values can then be utilised to determine the effective fractured porosity, which specifically focuses on the hydrodynamic behaviour of that particular segment of the rock formation and can be assessed alongside permeabilities in a consistent manner. By fitting power function curves on the standard poro-perm diagram, fractured intervals with similar fracture closure characteristics can be identified. Effective porosity and permeability values derived from the same dataset and related to the same rock body can also be used to calculate other important indices, such as the FZI mentioned above. The analysis of data from ten wells penetrating the fractured Mórágy granitoid in SW Hungary demonstrates that, based on the estimated hydrodynamic behaviour of individual intervals, hydraulic flow units with distinct characteristics can be distinguished in each well.

One of the main challenges in studying fractured reservoirs is extending porosity and permeability data to unexplored areas of the rock formation. Although a detailed discussion of this issue is beyond the scope of the present paper, having reliable fracture aperture coefficients for individual wells enables interpolation of aperture data between wells. Using actual hydraulic aperture data may provide a practical method for estimating viable porosity and permeability fields.