Fractal Characterization of Multiscale Fracture Network Distribution in Dolomites: Outcrop Analogue of Subsurface Reservoirs

: Fractured aquifers, especially dolomites, are important hydrocarbon reservoirs and sources of thermal and groundwater in many parts of the world, especially in the Alpine-Dinaric-Carpathian region. The most dominant porosity type is fracture porosity, which acts as the preferential ﬂuid pathway in the subsurface, thus strongly controlling ﬂuid ﬂow. Outcrops provide valuable information for the characterization of fracture networks. Dolomite rock properties and structural and diagenetic processes result in fractured systems that can be considered fractals. The fracture network was analyzed on 14 vertical outcrops in 35 digitized photographs. The values of the fractal dimensions varied slightly by the software and method used, but the trends were consistent, which conﬁrms that all methods are valid. Small values of fractal dimension indicate the dominance of a few small or large fractures, and high values of fractal dimension result from a combination of large numbers of small fractures accompanied by a few large fractures. The mean value of the fractal dimension for analyzed fracture networks was 1.648. The results indicate that the fracture network of the Upper Triassic dolomites can be approximated by fractal distribution and can be considered a natural fractal, and values can be extrapolated to higher and lower scales (1D and 3D).

Direct observation of the fault and fracture systems is often not possible due to the depths at which the rocks are located when they represent reservoirs of oil and gas, geothermal water, and even groundwater.In the petroleum industry, the only exact input data on fractures is one-dimensional (borehole data), but this is often insufficient to model the three-dimensional (3D) distribution of fractures [32].Also, borehole data refer mainly to small fractures representing a close area around the borehole [32,33].Discontinuity distribution and fault displacements (10 and more meters) are available after seismic data interpretation [32,34].There is a significant gap in knowledge of fault and fracture systems between the borehole and seismic data resolution.The characteristics of the fracture systems in this interval are the most important for characterizing the properties of the reservoirs or aquifers [32].To precisely model such fractures, a detailed survey of rock outcrops on the surface ("analog outcrop studies") can be of great help [7,[35][36][37][38][39].Rock outcrops are extremely reliable data sources on fracture systems because they contain continuous 2D to 3D data on their geometry and distribution, depending on some controlling factors such as the thickness of layers and the geological structure in which the observed outcrop is located [32,40].However, three-dimensional porosity data are rarely available, and their acquisition is often very expensive [41].Two-dimensional (outcrops, quarries, orthophotos, geophysical data, microscopic preparations, etc.) and one-dimensional data (borehole data, "scanlines", etc.) are much more common.
This research systematically investigated the fractal characterization of fracture system patterns of Upper Triassic dolomites in the Žumberak Mts. in northern Croatia.We digitized 35 outcrop images from 14 outcrops (Figure 1) to measure fractures in multiple scales.The box-counting method of fracture networks is the most popular fractal analysis method [66].These rocks are regionally the most widespread geological formation in the Alpine-Carpathian-Dinarides region, with huge significance for the industry since they represent a major geothermal aquifer.Furthermore, Upper Triassic dolomites are significant groundwater aquifers and hydrocarbon reservoirs that can be used for crushed stone [89].Due to intensive fracturing and/or recrystallization of these dolomites on the surface, there are only a few papers dealing with the porosity description and quantification, so this paper presents a significant contribution to the understanding of the porosity distribution and quantification of these rocks.

Geological Setting
Upper Triassic dolomites are one of the most dominant and regionally widespread geological formations in the Alpine-Carpathian-Dinarides region.Most of this Upper Triassic dolomite succession is presented by the Main Dolomite Formation (original name is Hauptdolomite [91].The most important characteristics of the Upper Triassic dolomites are huge lateral, widespread and large thicknesses from 300 m up to 3000 m [7,89,90,[92][93][94][95][96].These geological settings and significant porosity result in large groundwater and geothermal water aquifer potential. The research area is located in the Žumberak Mts, in the northwest part of Croatia, where the Upper Triassic dolomites are developed to a total thickness of 1590 m (Figure 1) [89,90].This area was selected because it is a nearly tectonically undisturbed area required for pore size analysis.Three formations were singled out: Slapnica (with members Vranjak and Drenovac), Main Dolomite (with member Kalje), and Posinak (Figure 1).Within the formations and members, cyclically alternating microfacies were determined: dolomicrite, fenestral dolomite, and stromatolitic dolomite [8,90].Sedimentary environments have been interpreted as a dynamic shallow-water system from subtidal to shallow intertidal [90].Tectonic settings are relatively simple, with few major faults, but rocks are not very tectonically disturbed (Figure 1).Layer inclination varies from horizontal up to 57 • in the vicinity of the faults.Most of the area is tectonically not disturbed, which was one of the main reasons for having preserved outcrops for fracture pattern analysis.

Theoretical Background
The concept of dimension in mathematics can be approached from different points of view, and there are different definitions of dimension [97].Euclidean and topological dimensions are finite and integer (1,2,3) and roughly define space.Many objects and phenomena in nature can be described with 0, 1, 2, or 3 dimensions (Figure 2A).However, Mandelbrot (1983) [98] noticed that many irregular objects and phenomena appear in nature that cannot be completely defined by topological dimensions.Examples are clouds that are not spheres, flashes of lightning that are not straight lines, mountain belts that are not cones, etc. [98].Their dimension is described by the Hausdorff (some authors also call it Hausdorff-Besicovich) dimension, which measures how much an object fills a metric space.To understand fractals and fractal dimensions, it is necessary to define concepts such as Hausdorff measure and dimension [99]: The δ-cover of a set F is a countable or finite union of sets {U i } with radii, 0 < |U i | ≤ δ that cover the set F. Let F ⊂ R n and let s ∈ [0, ∞ , for each δ > 0, we define: Therefore, the quantity H s δ (F) presents the minimal sum of the s-th powers of the diameters over all covers of F by sets of diameters not greater than δ.
As δ decreases, the set of permissible covers of F (in Equation ( 1)) is reduced [99], and hence, the value of the infimum H s δ (F) does not decrease and approaches the limit value as δ tends to zero [99] (Equation ( 2)): This limit exists for every subset F of the set R n , although the limit value is often 0 or ∞.H s (F) is called the s-dimensional Hausdorff measure of a set F. The Hausdorff measure is a measure that assigns some value from the interval [0, ∞] to each set in R n .A 0-dimensional Hausdorff measure of a set is the number of points in this set (if the set is finite) and ∞ if the set is infinite.The one-dimensional measure of a smooth curve or line in R n is equal to the length of that curve.This means the Hausdorff measure generalizes count, length, area, and volume but does not necessarily have to attain positive integer values [99].According to Equation (1), for an arbitrary set F ⊂ R n and if δ < 1, H s δ (F) does not increase with increasing s, and according to Equation (2), H s (F) will not increase either.Moreover, if t > s and {U i } is a δ-cover of the set F, then it holds: Calculating the infimum over all δ-coverings in Equation (3) results in Equation ( 4): By letting δ → 0 in Equation ( 4), it follows that if H s (F) < ∞, then H t (F) = 0 for every t > s.The graph shows that there is a critical value of s at which H s (F) "jumps" from ∞ to 0 (Figure 2B).This critical value is called the Hausdorff dimension of the set F, denoted dim H F, and is defined for each set F ⊂ R n : If the upper limit (supremum) of the empty set is stated value 0, then: If s = dim H F, then H s (F) can be 0, ∞, or it can satisfy: The Hausdorff dimension is often used to describe the fractal dimension (F dim ) of complex objects in nature.A point has Hausdorff dimension D = 0, a line D = 1, a surface D = 2, and a cube D = 3 (Figure 2A).Mandelbrot defined a fractal as a set, an object, or a phenomenon whose Hausdorff dimension (D) is greater than the topological one, where the dimension D is defined by [58,98]: where: N i -number of objects or fragments characterized by linear dimension r i ; Cproportionality constant; D-fractal (Hausdorff dimension) dimension that is calculated by (from Equation ( 8)): (9) where: ln is logarithm with base e, and log is logarithm with base 10.
Natural (statistical) fractals have the following properties [61]: (1) The parts of the object have the same structure as the object as a whole, except that they are slightly deformed in different scales (there are small fluctuations in the degree of fractality between scales)-a property of self-similarity.(2) Forms are often irregular and fragmented and remain so in all scales in which they exist.(3) They are created through an iterative process.(4) They have a fractal dimension.
The basic property of fractal distribution is scale independence.Mathematically, Equation ( 8) can be applied to an infinitely large-scale interval.However, natural fractals appear on a scale limited by lower and upper limits, which also applies to fractal distributions [100].Fragmentation should be considered to understand how to apply the fractal distribution to real data sets.Fragmentation is crucial to many geological processes and phenomena [61].Rocks can be fragmented in various ways-by tectonic processes that include faulting, the formation of fissures, and fracture systems through which the rock can be further fragmented by tectonic processes or weathering [61,100].The formation of fractures is a non-linear process, so even the simplest distribution of fractures and fracture systems is demanding and challenging to model.Fractal geometry allows quantifying geometric shapes that repeat over a certain range of scales [66].The fragmentation process (Figure 3A) and an example photograph of an Upper Triassic dolomite sample (Figure 3B) indicate that fractal geometry can quantify the scale range of the spatial distribution of fracture systems in dolomites.Also, the fractal description of fracture systems can represent a basis for analyzing fluid flow through rocks because the fluid flow that occurs through a fractal medium should also have fractal properties [43,59,66].The basic property of fractal distribution is scale independence.Mathematically, Equation ( 8) can be applied to an infinitely large-scale interval.However, natural fractals appear on a scale limited by lower and upper limits, which also applies to fractal distributions [100].Fragmentation should be considered to understand how to apply the fractal distribution to real data sets.Fragmentation is crucial to many geological processes and phenomena [61].Rocks can be fragmented in various ways-by tectonic processes that include faulting, the formation of fissures, and fracture systems through which the rock can be further fragmented by tectonic processes or weathering [61,100].The formation of fractures is a non-linear process, so even the simplest distribution of fractures and fracture systems is demanding and challenging to model.Fractal geometry allows quantifying geometric shapes that repeat over a certain range of scales [66].The fragmentation process (Figure 3A) and an example photograph of an Upper Triassic dolomite sample (Figure 3B) indicate that fractal geometry can quantify the scale range of the spatial distribution of fracture systems in dolomites.Also, the fractal description of fracture systems can represent a basis for analyzing fluid flow through rocks because the fluid flow that occurs through a fractal medium should also have fractal properties [43,59,66].

Materials and Methods
The Fdim of fracture patterns can be estimated from photographs of outcrops and samples (Figure 4) [43,52,100,[102][103][104].Image processing of the outcrop photographs can be a powerful tool for data acquisition, preparation, and analysis of fracture systems in rocks [103].As the photographs are characterized by two dimensions (length and width), the Fdim of the fracture systems in the photographs is 1 < D < 2. The applicability of Fdim in

Materials and Methods
The F dim of fracture patterns can be estimated from photographs of outcrops and samples (Figure 4) [43,52,100,[102][103][104].Image processing of the outcrop photographs can be a powerful tool for data acquisition, preparation, and analysis of fracture systems in rocks [103].As the photographs are characterized by two dimensions (length and width), the F dim of the fracture systems in the photographs is 1 < D < 2. The applicability of F dim in defining porosity features is that it theoretically allows the extrapolation of the estimated dimensions to a higher order (from 1D to 2D, or from 2D to 3D) [14,43,98,100,104] by a simple relationship:

Materials and Methods
The Fdim of fracture patterns can be estimated from photographs of outcrops and samples (Figure 4) [43,52,100,[102][103][104].Image processing of the outcrop photographs can be a powerful tool for data acquisition, preparation, and analysis of fracture systems in rocks [103].As the photographs are characterized by two dimensions (length and width), the Fdim of the fracture systems in the photographs is 1 < D < 2. The applicability of Fdim in defining porosity features is that it theoretically allows the extrapolation of the estimated dimensions to a higher order (from 1D to 2D, or from 2D to 3D) [14,43,98,100,104]   Equations ( 10) and ( 11) are valid because the intersection of a 3D fractal with a plane will result in a 2D fractal with the fractal dimension D = D3−D − 1.Both equations are valid only for isotropic and deterministic fractals [105].For real objects, in our case, dolomite rocks, the more applicable relation is [43]: Equations ( 10) and ( 11) are valid because the intersection of a 3D fractal with a plane will result in a 2D fractal with the fractal dimension D = D 3−D − 1.Both equations are valid only for isotropic and deterministic fractals [105].For real objects, in our case, dolomite rocks, the more applicable relation is [43]: Fractal analysis using computer programs is most often performed with the Boxcounting and radius mass methods.In the box-counting method, the photograph of the outcrop is covered with a square grid of different square dimensions.For each value of square size, the cells containing the pixels covering the object in the photo are counted.The F dim is then calculated by Equation (9).The computer program iteratively covers the given image with squares of different dimensions.The relationship between the size of the square and the number of squares that cover the object results in a log-log line whose slope represents the F dim of the object.The second most used method is the "Radius mass" method.With the "Radius mass" method, the algorithm arbitrarily determines specific points and draws circles of different radii that grow iteratively around them.The ratio of black (observed object) and white pixels (background) is counted at each iteration.Each step counts the total number of occupied pixels inside the circle.The paper used both methods for F dim calculation to compare results.
For the research, 14 outcrops in the Upper Triassic dolomites were photographed in detail, and 35 photographs were selected for further analysis (Figure 1, Table 1, Supplementary File S1).Numerous automatic algorithms were tested to digitize the fractures in the photographs.Still, the results were unsatisfactory due to large noise (fracture aperture and infill, vegetation, crushed rock fragments etc.) in the photographs, so software could not delineate fractures from the host rock.In other words, it was not possible to satisfactorily delineate fractures from the rock by automatic algorithms.Thus, all photographs were digitized manually (Figures 4 and 5).Photographs were adjusted in color, brightness, and hue to more easily distinguish fractures from the rock (Figure 5A).The outcrop was then delineated from the surrounding vegetation, and the region of interest (ROI) was defined.Fractures were then manually traced in the Adobe Illustrator 2021 software (Figure 5B).After all fractures were digitized, images were converted from color photographs into binary photographs (black and white) by ImageJ 1.53.csoftware to estimate the fractal dimensions of the fracture patterns (Figure 5C).The purpose of the procedure was to obtain binary photos that optimally showed the fractures in the rock, so that the software Fractal Fract.2023, 7, 676 9 of 18 operated only with white and black pixels.Processing is the most demanding process of photo preparation for fractal analysis.It includes all actions by which the image is prepared for conversion into the binary form: removal of vegetation, manual digitization of fractures, and conversion of the photograph into binary (Figure 5A-C).Binary photographs consisted only of black and white pixels; in our case, black represented fractures, and white represented rock (Figure 5C).After all fractures were digitized, images were converted from color photographs into binary photographs (black and white) by ImageJ 1.53.csoftware to estimate the fractal dimensions of the fracture patterns (Figure 5C).The purpose of the procedure was to obtain binary photos that optimally showed the fractures in the rock, so that the software operated only with white and black pixels.Processing is the most demanding process of photo preparation for fractal analysis.It includes all actions by which the image is prepared for conversion into the binary form: removal of vegetation, manual digitization of fractures, and conversion of the photograph into binary (Figure 5A-C).Binary photographs consisted only of black and white pixels; in our case, black represented fractures, and white represented rock (Figure 5C).In binary photos, if the fractures are represented by a value of 1, the rock is represented by 0. The vegetation is also represented by a value of 0, which reduces the Fdim of the fracture pattern in the photograph.To overcome the vegetation problem, it was necessary to create boundaries within the binary photo region of interest (ROI) within which the Fdim was evaluated (Figure 6).When the ROI was defined, the program calculated the Fdim only within the ROI, eliminating vegetation's influence on the calculation.
Several programs with several algorithms were used to compare the results and evaluate which program and algorithm gave the most reliable results.To estimate the Fdim of the fracture systems, ImageJ 1.53c®, FracLac (supplement for ImageJ) [106], and Fractalyse 2.4.® http://www.fractalyse.org,accessed 20 April 2021) were used (Supplementary File S2).Of the programs used, only FracLac could make an irregular ROI.An ROI could also be made using Fractalyse, but only in a square shape, which is often unsatisfactory (Figure 6).In binary photos, if the fractures are represented by a value of 1, the rock is represented by 0. The vegetation is also represented by a value of 0, which reduces the F dim of the fracture pattern in the photograph.To overcome the vegetation problem, it was necessary to create boundaries within the binary photo region of interest (ROI) within which the F dim was evaluated (Figure 6).When the ROI was defined, the program calculated the F dim only within the ROI, eliminating vegetation's influence on the calculation.For every outcrop photo, we calculated 2D porosity as the ratio between white and black pixels in ImageJ software.

Results
The calculated fractal dimensions differed between the software and methods used.Several programs with several algorithms were used to compare the results and evaluate which program and algorithm gave the most reliable results.To estimate the F dim of the fracture systems, ImageJ 1.53c®, FracLac (supplement for ImageJ) [106], and Fractalyse 2.4.® http://www.fractalyse.org,accessed 20 April 2021) were used (Supplementary File S2).Of the programs used, only FracLac could make an irregular ROI.An ROI could also be made using Fractalyse, but only in a square shape, which is often unsatisfactory (Figure 6).
For every outcrop photo, we calculated 2D porosity as the ratio between white and black pixels in ImageJ software.
(A) (B) Generally, the lowest values were related to fracture patterns of lower complexity, and higher values corresponded to higher complexity (Figure 8).Small values of Fdim indicated the dominance of a small number of large fractures or small number of small fractures (Figure 8A).Many smaller fractures, together with a small number of large ones, contributed to the complexity of the fracture system and, by that, to a higher value of Fdim (Figure 8B).The values for fractal dimensions obtained with all software and algorithms can be divided into three groups concerning the fracture system architecture: (1) Small Fdim values (Fractalysebox-count: from 1.465 to 1.525; Fractalysermass: from 1.592 to 1.650; and FracLac: from 1.492 to 1.600).Small values of fractal dimensions were obtained for outcrops with a relatively small number of fractures, i.e., less fragmented outcrops, or where large fractures dominated (Figure 8A).These were outcrops where not all sets of fractures had been developed or were not visible.In these parts of the outcrop, the fracture system's complexity was lower, resulting in lower Fdim values.(2) Average Fdim values (Fractalysebox-count: from 1.525 to 1.625; Fractalysermass: from 1.650 to 1.700; and FracLac: from 1.600 to 1.675).Most of the analyzed photos of outcrops belonged in this interval.This interval could be considered representative of the Upper Triassic dolomites of Žumberak.(3) Large Fdim values (Fractalysebox-count: from 1.625 to 1.670; Fractalysermass: from 1.700 to 1.802; and FracLac: from 1.675 to 1.725).The largest fractal dimensions were obtained for the more fractured parts of the outcrops, which resulted in a more complex fracture system (Figure 8B).Fractures filled large portions of the analyzed space Generally, the lowest values were related to fracture patterns of lower complexity, and higher values corresponded to higher complexity (Figure 8).Small values of F dim indicated the dominance of a small number of large fractures or small number of small fractures (Figure 8A).Many smaller fractures, together with a small number of large ones, contributed to the complexity of the fracture system and, by that, to a higher value of F dim (Figure 8B).The values for fractal dimensions obtained with all software and algorithms can be divided into three groups concerning the fracture system architecture: (1) Small F dim values (Fractalyse box-count : from 1.465 to 1.525; Fractalyse rmass : from 1.592 to 1.650; and FracLac: from 1.492 to 1.600).Small values of fractal dimensions were obtained for outcrops with a relatively small number of fractures, i.e., less fragmented outcrops, or where large fractures dominated (Figure 8A).These were outcrops where not all sets of fractures had been developed or were not visible.In these parts of the outcrop, the fracture system's complexity was lower, resulting in lower F dim values.(2) Average F dim values (Fractalyse box-count : from 1.525 to 1.625; Fractalyse rmass : from 1.650 to 1.700; and FracLac: from 1.600 to 1.675).Most of the analyzed photos of outcrops belonged in this interval.This interval could be considered representative of the Upper Triassic dolomites of Žumberak.(3) Large F dim values (Fractalyse box-count : from 1.625 to 1.670; Fractalyse rmass : from 1.700 to 1.802; and FracLac: from 1.675 to 1.725).The largest fractal dimensions were obtained for the more fractured parts of the outcrops, which resulted in a more complex fracture system (Figure 8B).Fractures filled large portions of the analyzed space irregularly, so the fracture network's complexity was high, resulting in high fractal dimension values.
Two-dimensional porosity (P 22 ) was calculated for 35 outcrop photographs and ranged from 2.84% up to 8.99%, with an average value of 5.7% (Figure 9).Since porosity is a threedimensional rock property, these values should be checked by some 3D method in the future to see how representative they are of the rock mass.The fractures were digitized as lines, and the aperture of the fractures was not considered, which could affect the porosity values.A correlation diagram between 2D porosity and F dim was prepared to analyze the dependence between these two parameters (Figure 10).There was no visible correlation between the two parameters, which indicated that more or less complex fracture systems can result in higher and lower porosities.Furthermore, for better fracture network description in 2D, it was necessary to analyze both parameters.
Fractal Fract.2023, 7, x FOR PEER REVIEW 13 of 20 irregularly, so the fracture network's complexity was high, resulting in high fractal dimension values.
Two-dimensional porosity (P22) was calculated for 35 outcrop photographs and ranged from 2.84% up to 8.99%, with an average value of 5.7% (Figure 9).Since porosity is a three-dimensional rock property, these values should be checked by some 3D method in the future to see how representative they are of the rock mass.The fractures were digitized as lines, and the aperture of the fractures was not considered, which could affect the porosity values.A correlation diagram between 2D porosity and Fdim was prepared to analyze the dependence between these two parameters (Figure 10).There was no visible correlation between the two parameters, which indicated that more or less complex fracture systems can result in higher and lower porosities.Furthermore, for better fracture network description in 2D, it was necessary to analyze both parameters.The fractal dimension of fracture systems is a measure of the complexity of the fracture system.It depends on the rock type and the tectonic history of the investigated area [43].Two-dimensional (2D) fractal dimensions close to 1 indicate a small number of frac- The fractal dimension of fracture systems is a measure of the complexity of the fracture system.It depends on the rock type and the tectonic history of the investigated area [43].Two-dimensional (2D) fractal dimensions close to 1 indicate a small number of fractures, i.e., a simple fracture architecture.Large values, close to 2, indicate an extremely complex fracture system: the rock is extremely fractured or fragmented by a complex network of fractures [101].The calculated values for the fracture systems in the photographs and the mean values of individual outcrops were mostly greater than 1.585 (Figures 7-10).
The mean value of F dim for all outcrops was 1.646.Dolomites are among the most fractured rocks in nature [6,43,107], and large values of the fractal dimension of fracture systems are expected of them.As the fractal dimension measures complexity, different fracture systems of different geometries but the same complexity will result in the same value of fractal dimension.The results indicated that the fracture system of the Upper Triassic dolomites could be approximated by fractal distribution and considered a natural fractal.Because of this, the spatial and statistical geometric properties of the fracture systems could be extrapolated to higher and lower scales by Equation (12).
The extrapolation of the fractal dimension from 2D to 3D by Equation ( 12) [43] resulted in F dim values for Upper Triassic dolomites, ranging from 2.53 to 2.86 (Table 2) depending on the software and the algorithm.Average values were between 2.62 (Fractalyse box counting) and 2.74 (Fractalyse Radius Mass) (Table 2).Generally, the results of this research were very consistent with results from [43], who analyzed similar rocks with similar methodology in Slovenia.

Discussion and Conclusions
The fractal dimension is a measure of the complexity of the analyzed pattern.Different fracture systems of different geometry but of the same complexity will result in the same fractal dimension.In this paper, we used manual fracture extraction from photographs.There are many semi-automatic or automatic algorithms published in the literature.Still, after trying a few of them, the results were not satisfactory due to several issues: (i) dolomites are highly fractured by fractures of multiple orders (from larger to smaller); (ii) fractures sometimes had an aperture that was filled with crushed dolomite or clay that was problematic for the algorithms to overcome; (iii) due to high fracturing, there were pieces of rock on the outcrop that were often considered as fracture traces, so the postprocessing of the results could take a long time.Considering everything, we found that only manual digitization (that could be biased) applied to this type of rock.
We analyzed the impact of fracture patterns in outcrop images on their fractal dimensions.The results can be summarized as follows: Fractal dimensions of fracture systems in outcrops can be subdivided into three categories: (1) Small 2D F dim values (Fractalyse box-count : from 1.465 to 1.525; Fractalyse rmass : from 1.592 to 1.650; and FracLac: from 1.492 to 1.600).(2) Average 2D F dim values (Fractalyse box-count : from 1.525 to 1.625; Fractalyse rmass : from 1.650 to 1.700; and FracLac: from 1.600 to 1.675).Most of the analyzed photos of outcrops were in this category.(3) Large 2D F dim values (Fractalyse box-count : from 1,625 to 1.670; Fractalyse rmass : from 1.700 to 1.802; and FracLac: from 1.675 to 1.725).
Simple Equation ( 12) could extrapolate from 2D F dim to 3D, and the results for the analyzed outcrops ranged from 2.53 to 2.86, with average values between 2.62 (Fractalyse box-count ) and 2.74 (Fractalyse Rmass ) (Table 2).
2D Porosity (P 22 ) for Upper Triassic dolomites was in the interval from 2.84% up to 8.99%, with an average value of 5.7%.
Dolomites are among the most highly fractured rocks [6,43,107], and high values of the fractal dimension of fracture systems are expected.As the fractal dimension is a measure of complexity, different fracture systems with different geometries but of the same complexity will result in the same fractal dimension.That is why the spatial and statistical geometric properties of the fractal fracture network can be extrapolated to higher and lower scales.Verbovšek, 2009 [43,108] analyzed the fractal dimensions of fracture systems on the Middle and Upper Triassic dolomites of Slovenia using different measurement methods.The values of the fractal dimensions for the "box-flex" and "box-rotate" methods were 1.78±0.02.The results of the "Box-counting-Full Method", which are comparable to the results in this research, were in the interval from 1.47 to 1.59 [43], which is, on average, somewhat lower than the results of this research.The porosity of the T 3 dolomite of the Sella Group (Dolomites, northern Italy) [7] can be compared with the results of this study.The fracture porosity of the Main Dolomites formation of the Sella group reached 2.4%.Such a "small" value of fracture porosity can be attributed to only two sets of fractures in the analyzed area [7].Published dolomite porosity values range from 1 to 27%, where porosity is not divided into fracture, sedimentary and diagenetic [6,7].This study makes a significant contribution to the analysis of fracture networks of dolomites (especially Upper Triassic dolomites) since there are only a few papers dealing with the subject due to excessive fragmentation and fracturing of dolomites, which make analysis almost impossible.
In carbonate reservoirs, the permeability is usually influenced by the fracture network [63].Depending on the fracture network characteristics, fractures can behave as fluid flow conduits, barriers, and combined conduit-barriers [63,[109][110][111].Detailed fracture description and parameter quantification then become crucial to analyze and predict the fluid flow behavior in subsurface aquifers or reservoirs [63,112,113].The results indicate that the fracture system of the Upper Triassic dolomites can be approximated by fractal distribution and can be considered a natural fractal.Because of this, the fracture system's spatial and statistical geometric properties can be extrapolated to higher and lower scales.Natural fracture networks in the subsurface are three-dimensional phenomena.With current technology, such as outcrop mapping, borehole logging, and 3D seismic surveys, performing a 3D fracture network survey is almost impossible.Therefore, investigating fractal characteristics of 3D fracture networks is challenging but very useful for fracture characterization and can be great input for DFN models that are used for 3D fluid flow modeling.

Fractal 20 Figure 2 .
Figure 2. (A) Fractal dimensions of Euclidean objects (point, line, area, and body with D = 0, 1, 2, and 3) and an example of a fractal fracture network with a fractal dimension of D = 1.76.(B) A graph showing the dependence of the Hausdorff measure ( ()) on the parameter s.The Hausdorff dimension is the value of s at which s falls from ∞ to 0 [99].

Figure 2 . 20 Figure 3 .
Figure 2. (A) Fractal dimensions of Euclidean objects (point, line, area, and body with D = 0, 1, 2, and 3) and an example of a fractal fracture network with a fractal dimension of D = 1.76.(B) A graph showing the dependence of the Hausdorff measure (H s (F)) on the parameter s.The Hausdorff dimension is the value of s at which s falls from ∞ to 0 [99].Fractal Fract.2023, 7, x FOR PEER REVIEW 7 of 20

Figure 3 .
Figure 3. (A) Theoretical model (three-dimensional fractal model of fragmentation) according to [101] and (B) fracture porosity in the Late Triassic dolomites in Žumberak Mts.
) whereby: D 3−D = fractal dimension of the object observed in 3D space; D 2−D = fractal dimension of the object observed in 2D space; D 1−D = fractal dimension of the object observed in 1D space.

Figure 3 .
Figure 3. (A) Theoretical model (three-dimensional fractal model of fragmentation) according to [101] and (B) fracture porosity in the Late Triassic dolomites in Žumberak Mts.
D3−D = fractal dimension of the object observed in 3D space; D2−D = fractal dimension of the object observed in 2D space; D1−D = fractal dimension of the object observed in 1D space.

Figure 4 .
Figure 4. Outcrop image processing and fractal dimension calculation flowchart.

Figure 4 .
Figure 4. Outcrop image processing and fractal dimension calculation flowchart.

Fractal 20 Figure 6 .
Figure 6.(A) Photograph of the outcrop with vegetation (FRA-9, IMG_PAN-51), (B) Example of ROI determination where vegetation occupies a significant area in the photograph.

Figure 6 .
Figure 6.(A) Photograph of the outcrop with vegetation (FRA-9, IMG_PAN-51), (B) Example of ROI determination where vegetation occupies a significant area in the photograph.

Figure 8 .
Figure 8.Comparison of fracture network with low fractal dimension (A) and high fractal dimension (B).

Figure 8 .
Figure 8.Comparison of fracture network with low fractal dimension (A) and high fractal dimension (B).

Figure 9 .
Figure 9. Two-dimensional porosity of the Upper Triassic dolomites obtained from outcrop photographs.Figure 9. Two-dimensional porosity of the Upper Triassic dolomites obtained from outcrop photographs.

Figure 9 . 20 Figure 10 .
Figure 9. Two-dimensional porosity of the Upper Triassic dolomites obtained from outcrop photographs.Figure 9. Two-dimensional porosity of the Upper Triassic dolomites obtained from outcrop photographs.Fractal Fract.2023, 7, x FOR PEER REVIEW 14 of 20

Figure 10 .
Figure 10.Relation of 2D porosity and fractal dimension of fracture networks in the outcrop photographs.

Table 1 .
Placemarks, photograph names (photographs are presented in Supplementary File S1), and basic geological information.

Table 2 .
(12)age, median, maximum, and minimum Fdim were calculated for individual images and averaged for outcrops.F dim ext is F dim extrapolated from 2D to 3D by Equation(12).