Modelling microfracture geometry to assess the function of a karst system ( Vízfő spring catchment area , Western Mecsek Mountains , Hungary )

One of the current objectives of karst research is to understand the spatial dynamics of karstic processes that increase ordecrease pore volume. Although we can construct 3D numerical models, it is a complex, multi-step process. These modellingapproaches combine a dissolution algorithm, a flow and/or transport model and an algorithm to reconstruct thespatial geometry of fracture networks. The paper focuses on the last task, and does not consider the first two problems. We applied the RepSim code, a DFN (discrete fracture network) type fracture geometry modelling software that usesfractal behavior of fracture patterns, for simulations. This method simulates fracture systems at a reservoir scale. The input parameters (length distribution, aperture, orientation and fractal dimension of the fracture midpoints) weredetermined using field measurements and evaluation of digitized images. The primary images were of outcrops fromthe surface and from two caves. The examined area is the karstic block of the Mecsek Mountains in SW Hungary, which is one of the most exploredkarstic regions in that country. There are seven small and one relatively large (Vizfő) catchment areas in the mountains. The large catchment was used as the study area, which lithologically consists of sandstones and limestones,both intensely fractured by subsequent tectonic events. The spatial distribution of cave entrances and dolinas is unevenacross the study area. This phenomenon has not yet been investigated. Here, a relationship is inferred between the original (prekarstic) microfracture network geometry and the spatial distributionof the aforementioned karst forms. The results show that the region can be divided into two zones that fracturedin distinctly different ways. Their fracture network communication features, porosity and permeability differ.Additionally, sub-regions could develop inside the catchment area where dissolution-cementation processes couldhave been differently effective, determining the spatial distribution of the dissolution governed karstic forms (e.g.caves, dolinas).

with the growing energy demand, the hydrocarbon industry is paying more attention to the exploration of fractured karstic reservoirs (LUN et al., 2010;FENG et al., 2013).
Karstic evolution is a complex process, and the development of these areas is influenced by many factors in both space and time, which may result in an extremely heterogeneous structure and morphology.The quantity and quality of precipitation, the ever-changing surface and the subsurface geomorphology and biological processes are equally important factors.Lithological and structural parameters are important as they are more constant in space and time those biological or environmental factors.In numerous limestone bodies where the syngenetic primary porosity has drastically decreased during diagenesis, the main flow paths are fractured structures that developed through different deformation events.They are linked to other "inception" horizons, whose size, form and permeability can greatly change during karst formation.According to PALMeR (1991), 57% of karstic conduits are related to bedding planes, 42% to fractures and 1% to the rock matrix.
One of the main objectives of karst research is to understand the spatial dynamics of karstic processes that increase or decrease pore volume.In the early 1990s, karst system modelling was mainly used to reconstruct individual fracture developments (e.g., PALMEr, 1991;DrEybroDT & GA-broVšEK, 2000b;bAUEr et al., 2000).From the mid-1990s, technology enabled 2D models to also be used (KAuF MANN & brAUN, 1999, 2000;bAUEr et al., 2002bAUEr et al., , 2003;;KAUFMANN, 2003a,b;DrEybroDT et al., 2005).We can now construct 3D numeric models, even though this is a complex, multi-step process.Three-dimensional numerical models are typically based on three key aspects: a dissolution algorithm, a flow and/or transport model, and an algorithm to reconstruct the spatial geometry of fracture networks.Our paper focuses on the last task.
The aforementioned studies, depending on their purpose, were interested in reconstructing initial fracture geometry (e.g., JACqueT et al., 2004, KAUFMAN et al., 2010;borGHI et al., 2012;PArDo-IqúzUIzA et al., 2012).Several papers discuss some regular structure or formula to which the spatial net is applied.The development of computer technology has made it possible to realistically reconstruct more complex geometries, such as karstic networks.The fundamental precondition of a near-reality reconstruction is that the input parameters of the model should be based on measurable factors from the known area.Karstic objects are fractal-like objects (e.g., length distribution of caves, (CUrL, 1986)), so their geometric parameters given in a particular size range can be scaled to both smaller and larger size ranges.It is also notable that, because the quantity and quality of the perceptible factors are variable within a given range, very small and very large forms may not necessarily fit spatial reconstructions.It is best to sample the range of cores, hand specimens and cave passages, and the latter was where PArDo-IqúzUIzA et al. (2012) began.because more geometric information can be obtained by analyzing microfracture networks in less developed or speleologically less explored karsts (e.g., karstic hydrocarbon reservoirs), we surveyed the parameterizability of these fractures.
Individual fractures that co-occur, such as karstic forms, can appear in a wide range of scales from submicroscopic to hundreds of metres in diameter (ALLègRe et al., 1982;oUILLoN et al., 1996;TUrCoTTE, 1997) and their features can be compared across different ranges.The structural data of the fractures and microfractures obtained in surface explorations, caves and boreholes can provide important information about large-scale structures, faults and fault-zones.The qualitative features and the values of the measurable parameters (size, spatial density) depend on the petrographic features of the given deformed rock type, e.g., its mineral composition, grain size distribution, structure and former orientation.The physical circumstances in which the deformation takes place also play a vital role; the primary significance of the stress field is substantially modified by the lithostatic pressure, temperature and the rock's fluid content.The resulting fracture system, as a geometric object consisting of individual fractures, is created as the cumulative result of all these forces.Although under the Coulomb criterion it is possible to prove the genetic connection between the stress field and the geometry of the fracture system, the former still cannot be deduced from the latter (CAMPOS et al., 2005).Therefore, the simulation algorithm (dealt with geometry), unlike several others (e.g., oLSEN, 1993;rENSHAW & PoLLArD, 1994;rILEy, 2004), does not examine the evolutionary process of the fracture network in the given stress field, but aims to reconstruct the complex geometry.
The karstic block of the Mecsek Mountains in SW Hungary has been investigated for several decades.The purpose of our paper is to supplement the experimental (dye tracing, well data, speleological observations, etc.) modelling results of the fracture system to help understand how the main karst system of the mountains (Vízfő) functions.

GeOlOGIcAl bAcKGrOUND
The most important rock types for karstification in the western Mecsek Mountains are the diverse Triassic formations, particularly the Anisian Lapis Limestone.The oldest Triassic formation in the area is the Jakabhegy Sandstone Formation (bArAbáS & bArAbáSNé, 2005), which is composed of variations of quartz-cemented, cross laminated polymicritic sandstones and conglomerates, followed by the Hetvehely Dolomite Formation, with a total thickness of not more than approximately 40 metres (HAAS et. al., 2002) .Next is the Lapis limestone, a highly fractured rock, with relatively high clay and organic content; this formation is thinly layered, sometimes bedded and formed in an innerramp facies.because of the syn-sedimentary deformation events (syn-sedimentary faults) and post-sedimentary faults, the original stratification is highly disturbed and can only be traced in blocks of some tens of metres.
because of Alpine orogenic movements during the Cretaceous, the whole sequence is strongly folded (bENKoVITS et al., 1997;CSoNToS et al., 2002).The study area is located on the northern limb of the main fold (Fig. 1).The Mecsek Mts. were raised and thrust on the northern foreground and partially folded the young sediments (SEbE et al., 2008).
Athough the present structure of the mountains was determined by these Mesozoic movements, the karstic processes were more strongly influenced by younger tectonic events.
The NE-SW structural lines, where fissure caves and karstic cavities formed, opened in the Pleistocene in several stages (SzAbó, 1955).After this tectonic event, N-S and NW-SE structural lines formed.In the early Holocene, a general rise in the region occurred, which resulted in a significant increase in erosion (SzAbó, 1955).Simultaneously, the karstic water level began to heavily fluctuate, which influenced the valley structure of the western Mecsek Mts. and strongly affected the formation of its caves (SzAbó, 1955;LoVáSz, 1971;SEbE et al., 2008).The present direction of the largest main stress is Ne-SW, which opens parallel fractures and closes perpendicular ones (bADA et al., 2007) and has a strong effect on the relationship of the fracture system.The karst-forming layers now typically dip 20° to the north and are heavily folded.No thorough survey has been prepared on the hydrogeological impacts of the structural elements.
The neotectonic and geomorphological evolution resulted in eight catchment areas in the central karstic region in the western Mecsek Mts.(róNAKI, 1970) (Fig. 1), the largest of which is the catchment area of the Vízfő spring.Approximately 30% of the 16 km 2 catchment area consists of non-karstic formations; the Jakabhegy sandstone is the highest relief in the south, while the Lapis limestone domi-nates in the north.Our research area covers approximately 60% of the whole catchment, covering the territory between the Szuadó-valley, the Körtvélyesi-valley and the Vízfő cave (Fig. 2).Data from the most important caves are summarized in Table 1.
According to SzAbó (1955), the karst phenomena in this catchment area formed in several individual stages and the valley formation and the subsurface karsts started to develop simultaneously.There are 73 discovered cave entrances in the 16 km 2 drainage area.After several decades of speleological research, it is now clear that most of the caves and cave parts (passages, halls) are oriented to the northeast (Fig. 3), and numerous caves have siphon endpoints.The  spatial distribution of these caves, however, is uneven, similar to that of dolines.LIPMANN et al. (2008) found that dolines are concentrated in the central part of the area, defining NE-SW oriented chains for structural reasons.The shapes of individual dolines was not determined by structural effects, but rather by micromorphological and microclimatic circumstances.The average annual precipitation is approximately 600-650 mm, (higher than the Hungarian average), which is the result of the strong Mediterranean effect.The least precipitation occurs in January, and the most occurs in June, although there is a second peak in the autumn.
There are two permanent surface water courses in the area: the orfű-creek in the Szuadó-valley and the Körtvélyes-creek in the Körtvélyesi-valley.Their average rate of flow is 52 m 3 /day and 17 m 3 /day, respectively.The water from both creeks is directed to the Vízfő spring through sinkholes, but our experience shows that a large amount of water leaks underground through the fracture network from the creek beds.Little information has been found on the path of the water through the sinkholes.The test results are shown in Table 2.The average flow rate of the Vízfő spring is 420 m 3 /day, but can exceed 100,000 m 3 /day during a flood.(modified KONRÁD et al., 1983).

MethODs
Individual fractures can be interpreted as 2D surfaces of finite extension, usually bent in multiple and very complex ways, but can mostly be approached as planar (CHILES & DE MArSILy, 1993).In our case, we followed the circle representation in both the parameterization of fractures and further simulation.The geometric parameters clearly describe individual fractures as the spatial position of the circle's center, the radius and orientation (dip, strike).For a fracture system, this can be interpreted as a function of the spatial density of the center points and a distribution function typical of the radius and dip-strike value-pairs.The hydraulic description of the fracture system needs the positive volume of the individual fractures; therefore, circle plates without thickness were replaced by flat disks with apertures (parallel plate model, WITHErSPooN et al., 1980).Discrete fracture network modelling requires characterizing the geometry of individual microfractures, which was completed using the methods described below.Abbreviations used to describe the fracture network geometry are summarized in Table 3.

Field work
We first localized all outcrops on the surface and subsurface, recording their gPS coordinates.The orientation (α, β) of each individual fracture was measured using a compass after the outcrops had been cleaned.Next, 5 MP photos were collected from representative sections of the outcrops.In 9 cases, hand specimens were also collected to investigate their microfracture systems under laboratory conditions.using the photos, thousands of individual fractures were digitized from hand specimens using ArcView software.Digitized images were used to detect the length, distribution, aperture and the relative position of the fractures Neither the shortest nor the longest fracture sets fit on the expected line.both biases come from representative problems at extreme sizes, so they were left out of the subsequent analysis.Outlier data were selected using the grubbs test.
Measuring the fracture aperture is an uncertain task.Although there are accurate and accepted methods in the literature (e.g., MicroCT, MrI), these techniques require special laboratory equipment.Applying these approaches is impossible in the field.In many cases, the distance between the opposite fracture walls could not be measured correctly because of the irregular and rough walls.The aperture of the individual fractures was measured at least three times at three different points (typically four or five).To evaluate the data, the average values of the measured data were used.
Although there were many uncertainties, previous investigations suggest that the aperture follows a power law distribution similar to the length (DE DrEUzy et al., 2002;ORTegA et al., 2006).In addition, a linear correlation exists between the length and the aperture of the individual fractures (bArToN & LArSEN, 1985;VErMILyE & SCHoLz, 1995;GUDMUNDSSoN, 2000;GUDMUNDS-SON et al., 2001).The aperture of a single fracture can be described by the following equation: The two parameters (A, B) can be determined by computing the linear regression between the measured length and the aperture data.
The spatial distribution of fractures was analyzed using the fractal dimension of the midpoints of the digitized fractures.There are numerous methods to calculate the fractal di-

Mapping
Fracture parameter maps were generated to spatially extend the data measured in discrete points to the whole study area.All maps were computed using Surfer 8.0 using the linear kriging approach because there were not enough points to thoroughly fill variography calculations.The grid cells were 40x40 m.

Modelling
The RepSim code was used many times for 3D fracture network simulations (M.TóTH, T. et al., 2004, VASS, I. et al., 2009).repSim code has been applied for building the DFN (Discrete Fracture Network) model because it follows the mension, such as Mass radius, Cumulative Intersections, Vectorized Intersections, Convex Hull Intersections, Convex Hull, etc.We applied the box Counting algorithm (MANDeL-broT, 1983;MANDELbroT, 1985) to characterize the fracture network following the suggestions of bArToN & LARSeN (1985) and bArToN (1995).The benoit 1.0 software was used, which demands 1 bit (black and white) images of the fracture midpoints.Two breakpoints typically appear on each point set generated using the box Counting algorithm.Points below the first one belong to very small box sizes, while those above the second one come from the largest boxes.The two point classes have no geological meaning and are only derived from the applied approach.A regression line was fit between the two breakpoints in each case.
box-counting algorithm to locate fracture centres and uses the measured D, e and orientation data for simulation.This is a DFN modelling software that uses the fractal geometry behaviour of fracture patterns and is based on the previously mentioned parameters.Although many papers suggest that the assumption of circular shaped fractures is more likely to validate in massive rocks, there is absolutely no data to use a polygonal or ellipse approach for the individual fractures.Although the simulated system consists of a myriad of discrete flat cylinders, hydrodynamic parameters (porosity, intrinsic permeability tensor) are calculated for a given grid net for hydrodynamic modelling.The size of the grids is defined following the rEV (representative elementary volume) concept of bEAr (1972); the volume, above what the calculated porosity does not change remarkably.Following the approach of M. TóTH & VASS (2011), rEV is defined at the grid size where the variation coefficient of porosity becomes stable.It is important to note that porosity, permeability and ReV calculations relate only the fracture system and the model ignores the effect of the bedding planes.

DAtA
Measured fracture length data follow the expected power law distribution in each studied locality, so the E and F parameters needed for modelling were derived.Similarly, a box counting calculation was successful in each case, and the D parameter, representing the spatial distribution of fracture midpoints, was computed in all 11 study points.The results from the parameter calculations are in Table 4, and the evaluated data of a typical sample are in Figure 4.The results from the aperture measurements show some deviation from the assumed linear relationship, likely because of dis-solution for a few microfractures.Most fractures lie along the trend line (Figs.4E,F) if the strongly dissolved fractures were ignored during the linear regression calculations (c.f.Equation 2).Parameter maps calculated using normal krigings are illustrated in Figure 5.

simulated fracture network
At the beginning of modelling, we defined the primary cell edge as 10 m, which was halved during the iteration process until we reached a minimal cell length of 1 m.This method required four iteration steps, during which we generated approximately 113000 single fractures.For each midpoint, a sin-  gle fracture was chosen from the 1-50 m length interval using the calculated length distribution.Altogether, 10 independent fracture network models were computed and evaluated.
In each case, the generated model forecasted two fields of dramatically different behaviours: a highly fractured zone, which has at least one well-commutating fracture group (Type 1) and a zone, which has no well-communicating fracture cluster at all (Type 2) (Fig. 6).The connectivity of a fracture group can be defined by the proportion of the interconnected fractures.It essentially depends on the D, E and orientation values of the fracture network (M.TóTH & VASS, 2011).In this paper a fracture group was considered well-communicating if it contains more than 5% of all simulated fractures from the modelled volume.
Two large, independent fracture groups were rendered probable in the area.One on the southeastern part contained approximately 10% of all of the fractures (Type 1A), while another communicating system in the north contained approximately 5% (Type 1b).both zones had relatively low D values (~1.26) and relatively low E values (<1.00), suggesting that the fractures were relatively long.In the Type 2 zone, geographically surrounded by Type 1A and Type 1b areas, even the largest fracture networks consist of no more than approximately 0.1% of all fractures.This region is characterized by high D values (~1.6) and high E values, which suggests that, even if a dense fracture system was typical here because the network is dominated by short fractures, the number of communicating subsystems is subordinate (Fig. 7).
To test how well-connected the Type 1 systems are, each fracture with a length below a certain threshold was diminished to model the possible role of vein cementation.As the simulation of cementation goes on, the size of the communicating clusters decreases evenly.even if fractures with apertures of 120 mm (L = 20 m) are closed, a remarkably large communicating fracture group can be rendered probable, suggesting a mature 3D fracture system is present both in the northern and southern part of the study area.
The variation coefficient calculated from the average and standard deviation data of porosity for different grid sizes stabilizes at approximately 0.6, so we defined the rEV typical of Type 1 zones in the corresponding 50 m.The fractured porosity for fractures of at least 1 m (aperture 8 mm) for this rEV is 6.8%, which is considered good porosity for Triassic reservoirs.If we also include other possible porosity types (solution, along bedding, etc. (CHoqUETTE & PrAy, 1970), however, this area seems to have outstanding porosity.The permeability tensor calculated for the representative volume shows significant NE-SW anisotropy of 1.37 with an average value of 8*10 -11 m 2 .This direction does not depend on the orientation of the particular fractures but is closely connected to the anisotropy of the fracture density (Fig. 6).
In Type 2 zones, the biggest connected fracture group in the models only generates 0.01% of the fractures, while plenty of small fracture groups are evenly distributed.ReV is defined as 70 m for this area.The fractured porosity calculated for a volume above rEV is < 0.001%, which means that this zone has neither storage nor conducting capacities.That is why we believe that the slight Ne-SW permeability anisotropy typical of this area has no practical geological effect.The calculated average fractured permeability of 6*10 -13 m 2 must be even smaller because the fractures do not define a connected network.
According to the previous statements, one can see that the fracture system of the examined area and its hydrodynamic features (porosity, permeability) are not uniform.The northern and southern parts of the area have outstanding effective fracture porosity, while in the middle, it is practically negligible.
Considering the hydrodynamics and karst formation of the area however, it is extremely important to know the joint features of the three different zones.The rEV calculated for the whole area is approximately 90 m, which is rather high compared to other case studies (e.g., WU & KULATILAKE, 2012).This increase in size is because linking the two types of areas results in a substantial growth in the standard deviation of porosity.The fractured porosity above this ReV is 2%, which is an average effective porosity value for the  whole of the studied karstic reservoir.using other methods, VASS (2012) obtained similar results for the catchment of the nearby Tettye spring, where the lithological and structural environment is the same.

Geological and morphological consequences
Although the dominating Lapis limestone has different varieties, there is no evidence for different dissolution abilities.Additionally, current and palaeo-climatic parameters do not differ either; there are no places with local maximum and minimum values of precipitation.There are no substantial differences in the elevation of the study area, so we can assume that the main factors influencing epigenic karstic formation (soluble rock, precipitation of appropriate quality and quantity) affected the whole area to the same degree.There can be a crucial difference in the state of deformation of the rocks, which can affect the present surface and subsurface morphology as well.The modelling results suggest that the spatial distribution of the rate of brittle deformation is extremely heterogeneous in the area.
The Ne-SW orientation manifested in fracture density maps is also shown by structural geological surveys performed on a larger scale in the wider area of the Mecsek Mts.(CSoNToS et al., 2002;LIPMANN et al., 2008;KoNráD & SEbE, 2010).This renders the genetic link of micro-and macro-structures probable, so our observations are coherent with the statement that fracture patterns follow a scale-invariant pattern and can be modelled by fractal geometric approaches (yIELDING, 1992;KorVIN, 1992;NIETo-SA-MANIEGo et al., 2005;M. TóTH, 2010).
Although this research did not target the analysis of the fracture-system genetics, the spatial distribution of the different fracture systems is still closely related to well-known tectonic events.We found that when the short fractures in the model were blocked (e.g., mimicking cementation), the size of the fracture system in the southern part of the area (system Type 1A) decreased from south to north.Our measurements showed that the basic difference between the northern and southern regions of the study area, compared to the central region, was the length distribution of the fractures.One reason for such a crucial spatial variability of fracture lengths may be the complex tectonic evolution of the region.The deformation response to the given stress field systematically changed when moving away from the anticline axis of the western Mecsek Mts. (Fig. 1).Relatively few, but larger fractures formed near the axis, while a denser fracture system of smaller fractures formed (Type 1A) further away.The fracture system in the south seems to be fold related, but the northern region is related to the fault system that created the orfű basin (Type 1b) (Fig. 8).
GAbroVSEK et al. ( 2004) and SINGUrINDy & bErKoVITz (2005) observed that the net efficiency of karstic processes is significantly higher in well-fractured zones where the permeability and flow velocity are much higher inside a given rock mass.In these regional processes doline formation, karren formation, or even cave formation may be enhanced.In the less fractured zones, fractures are cemented and karsts typically do not form.In intensely fractured zones, a positive feedback process occurs.During this process, the more opened fractures have larger soluble surfaces, which enhance solution processes, so the karst formation process speeds up (JAKUCS, 1971;M. TóTH, 2003;FILIPPoNI et al., 2009).Structural control like this is wellknown in several karstic areas (bILLI et al., 2007).
The results suggest that the relationship between different karstic forms and the fracture systems should be independently studied for the three diverse fracture patterns.type 1 These areas can be characterized by a communicating fracture system with high permeability and porosity values.We can expect significant solution forms in this type of area.There is no difference in the karstic phenomena between Type 1A and Type 1b and all known large caves are all Type 1.These caves (e.g., Trió cave, Szuadó cave, Gilisztás cave, Vízfő cave, Fig. 2) formed through solution along fractures and erosion.Tectonic preformation substantially influenced and is influencing cave formation today.Although the regional flow gradient follows to the north because of the position of the anticline, all of the aforementioned caves are oriented Ne-SW, corresponding to the orientation of the permeability anisotropy previously calculated.Dolines are relatively rare in these two zones (LIPMANN et al., 2008).type 2 SINGUrINDy & bErKoVITz (2005) discuss that in areas of low porosity and permeability, cementation processes dominate.In the central zone of the study area, this statement is supported by the lack of significant (longer than 50 m) caves.The orientations of all known small cave-entrances and immature vertical caves do not follow the permeability anisotropy.The very low rate of underground karstic processes is because the zone has little fractured porosity and no ability to conduct water.
In this zone, the number of dolines is remarkably high (LIPMANN et al., 2008), but NE-SW orientations are atypical.As these dolines are neo-formations (LIPMANN et al., 2008), we think that the geometrical features of micro-fractures must have determined the initial stage of doline and doline-chain formation according to SINGUrINDy & bErKoVITz (2005).The relatively high number of individual fractures (D=~1.6)can be the main reason for the high doline number.Although the there are many dolines in this area, these are very immature probably as a result of the low permeability of the fracture system, which is coherent with the conclusion of SINGUrINDy & bErKoVITz (2005).Additionally, a significant number of the small caves in the area are linked to dolines from a solution origin.CzIGáNy and LoVáSz (2006) emphasized the definite NE-SW orientation of the doline chains without any orientation of the particular dolines.This agrees with the orientation of the large structural elements throughout the study area, while the microfracture networks are immature.

theoretical model
The two large fracture groups (Type 1) found in the southwestern and northern part of the study area can be identified with the spatial location of the most important caves.The model predicts the largest interconnected fracture system near the Szuadó valley caves, while the second largest interconneccted zone is predicted to be near the Vízfő cave (Fig. 9).In the case of the Szuadó valley sinks, the fracture model assumes a connected fracture system exists, which suggests that the Szuadó, Gilisztás and Trió caves are directly connected.This does not necessarily mean the connection can be passed through by people.Radon gas exploration was used for the Szuadó valley sinks (KoLTAI et al., 2010) to explore the connection.The series of measurements was not successful, as it could not prove or exclude connections.This issue still requires more research as earlier explorations did not target the relationship of caves to each other.There is a hydrological connection between the sinks and the spring areas (GILA, 2000), although no passages have been found between the two highly fractured regions.Only a certain amount of water leaking into the karst moves along the fracture planes.Without these planes, a large volu me of water may move along bedding planes (KAuFMANN, 2003;KAUFMANN & roMANoV, 2008).A typical example is the initial section controlled by planes in the Szuadó valley sink caves (e.g., Trió cave), which can carry as much as several hundreds of litres of water per minute (unpublished).In these caves, these sections are typically oriented north -south, so the strongly northerly average plane inclination on the northern wing of the anticline contributes to the northeast run of the caves.We observed that, after a short initial section of the sinks, the flow tracks connect to mature fractures.The flow rate does not grow along the planes during floods, but even small fractures have significant flow rates.Although the karst and spelaeology literature provides information on a large number of caves with important galleries following the bedding planes, the cave evolution in the Mecsek Mts.depends mostly on the fracture network and fault system.The bedding is disrupted by syn-and post-sedimentary faults, fractures and folds, so the lateral spreading is smaller than a few tens of metres.
Although there are no flow-place measurements (in fractures, faults or on bedding planes) of karstic water, passage features in the Vízfő spring between the sinks and the spring were examined by dye tracing tests (GILA, 2000;SzŐKE & orSzáG, 2006;róNAKI, 2007).Although these surveys have shown a short time before appearance of the dye and thus a strong hydrological connection, because of the very high deviation of the first arrivals, these passage feature estimations are extremely uncertain (Table 1).Despite the rapid appearance, each dye tracing test experiment showed some slowing features in the hydrological system.GILA (2000) and SzŐKE and orSzáG (2006) assumed that slowing by siphons or storage lakes occurs because the transition time is not stable; it can be short or long.The previously discussed model that predicts the coexistence of different fractured zones in the area implies that the slowing factor between the sinks and the springs should not be a series of siphons or storage lakes, but that the poorly connected fracture system is responsible for this phenomenon.The slowing rate of the middle zone must depend on the quantity of the received water.The extremely rainy period in 2010 is a good example of this.In 2010, there was a fluctuation of more than 50 m in the water level of the Szuadó valley caves, while in the other part of the karst, around the spring, the water level did not even reach 5 m (unpublished).Using the current fracture model, in the Type 2 zone, the fracture network is unconnected and the bedding planes alone are responsible for water flow.As soon as they became saturated, further tapping was impossible along bedding planes where there is low permeability.The middle zone then slows and diverts the flow from a northern to a northeastern direction.The fact that the karstic water arriving at the southern border of the non-communicating zone has to turn to the northeast is a very important element in the hydrogeological system and the development of the cave orientations.
To sum up, our model suggests that in the catchment of the Vízfő cave between the Szuadó valley caves and the Vízfő spring, a zone should exist where cave maturity is significantly underdeveloped compared to its vicinity because of tectonic preformation.This zone plays an important role in the hydrology of the area as it greatly decreases the permeability between the sinks and the spring, and practically diverts flowing waters to a northeasterly direction.

cONclUsIONs
• Through the study and modelling of the microfracture network in the examined area, the vicinity previously considered homogeneous can be divided into three distinctly different zones which differ in fracture network connectivity, porosity, permeability and type of karst development.• In the study area, the rEV calculated on fractured porosity is approximately 90 m, but for the well-fractured sub-areas, this volume is much smaller (50 m).• There are direct fracture geometric connections between the sinks of the Szuadó valley (Szuadó, Gilisztás and Trió caves), but they are not necessarily passable by people.• We have increased the knowledge of the karstic development in the area with new elements by studying the microfracture network and applying the previously discussed simulation method.This method makes it possible to simulate fractures on a reservoir scale using point-type data (even drills).This can be extremely important for exploration of hydrocarbon and geothermal reservoirs.

Figure 1 :
Figure 1: Location of the study area and the karst aquifers of the Western Mecsek Mts.(modified after RÓNAKI, 1971).

Figure 3 :D
Figure 3: Orientation of the largest cave passages in the study area (N=55).
based on previous studies, the distribution of fracture lengths (measured on 2D surfaces) follows a power law distribution, which can be approximated by a straight line on a log(N(l)) -log(l) plot(e.g., yIELDING et al., 1992; NIETo-SA- MANIEGo et al., 2005).Calculating the frequency distribution is the first step for the precise definition of the power law functions.because the number of cases on the histogram can have a great influence on its shape, this number was defined as k = 2 * INT(log 2 (max(l))

Figure 4 :
Figure 4: Processing methods of a typical sample.A, B) Result of the box-counting analysis of fracture midpoints.C, D) Frequency distribution of the fracture length values and the corresponding linear regression function.E, F) Relationship between the aperture and length values, and the same relationship without the extremely dissolved fractures (applied software: SPSS 20).

Figure 5 :
Figure 5: Spatial distributions of the measured parameters.a) Map of the D parameter.b) Map of the E parameter.c) Map of the aperture.White circle = sampling locations; black circle = cave entrances.

Figure 6 :
Figure 6: Calculated communicating fracture groups, and the horizontal segments of the permeability anisotropy ellipsoids in the three differently fractured domains.

Figure 7 :
Figure 7: Typical horizontal sections of the different types of the simulated fracture networks (not to scale).

Figure 8 :
Figure 8: Cross-section along the A-A' line (see Fig. 2) with the conceptual fracture systems.

Figure 9 :
Figure 9: Effect of the fractures and bedding on the hydraulic behaviour of study area.The main conductive zones are fault related fractures (yellow), and according to the DFN m del, there is an unfractured, bedding only, hydrological connection between the inflow and outflow regions.

table 1 :
Topographical data of the main caves in the study area.

table 2 :
Results of the tracing works in the Szuadó-valley.
Asl: Above Baltic Sea Precipitation: Difference of the previous three mounth precipitatin sum and the long term average precipitation az the same three mounthN.D: No Data available

table 4 :
Measured geometrical parameters of the fracture network.