Characterization of karstic aquifer complexity using fractal dimensions

The present study analyzes fractal dimensions for the daily discharge data series of 12 karstic springs registered over two decades in Northeast Hungary. Fluctuation in the observed data is frequent and irregular, producing rough time series. The level of roughness is measured by the fractal dimension defined in different ways and corresponds to the intensity of fluctuation. That, in turn, results from the structure of the karstic aquifer, its conduits’ geometry, and the water migration in them. In the given case of springs, p-variogram based fractal dimensions reflect the karstification level primarily. On the other hand, box-count and information dimensions are associated with mixing karstic and hydrothermal components when the latter is present. Therefore, the analysis of fractal dimensions of spring discharges may provide a way to obtain information on the complexity of the hidden subsurface conduits and the water flows in them in an exploratory and comprehensive way.


Introduction
Karstic catchments are characterized by spatial and physical heterogeneity, very complex fractured and porous structures, and water migration in these media appears to be highly non-linear. Therefore, linear models are unable to reconstruct the temporal variability of karstic spring discharge with sufficient accuracy. For a deeper understanding of the processes, other methods must be used. Several studies (Tessier et al. 1993(Tessier et al. , 1996Pandey et al. 1998;De Lima and Grasman 1999) demonstrated that a strong feature of rainfall and river discharge time series is their multifractal nature. Multifractal property characterizes flows through porous media (Veneziano and Essiam 2003), rock roughness, cave labyrinths (Andreychouk et al. 2013), and many other geological objects or phenomena.
To our knowledge, the multifractal properties of the discharge time series of karstic springs were first recognized in Labat et al. (2002), and subsequently further papers (Márkus and Kovács 2002;Márkus 2003;Majone et al. 2004;Livina et al. 2007;Labat et al. 2012) dealt with the fractal nature of karstic springs. In particular, in Márkus and Kovács (2002), two of the current authors reported on the multifractal character of some springs included in the present study. The word "fractal" originates in Mandelbrot's works, e.g., Mandelbrot (1977). He intended to describe the rough, broken, and irregular character of certain objects by this notion. It is characteristic of fractals that roughness is present at all scales. Fractals fill the Euclidean topological space-the one they are embedded into-more thickly than an object of lower topological (integer) dimension. However, they still do not fill out completely any open set of that Euclidean space. Hence they have a non-integer dimension, the so-called fractal dimension. The scaling properties of spring discharges may shed light on important hydrogeological characteristics of the corresponding aquifers. A comprehensive qualifier of scaling is the fractal dimension.
The present study considers the fractal character of the observed daily discharge time series of 12 karstic springs located in Northeast Hungary. The use of various FD estimation methods emphasizes various aspects of roughness, and consequently, significantly different FD values can be obtained.
The present paper analyses the box-count estimator, the power variogram (or pvariogram for short), i.e., the variogram, madogram, rodogram, and Increment1 based estimators (cf. Gneiting and Schlather 2004;Gneiting et al. 2010) for estimating and comparing the FDs of discharge time series.
The Rényi or generalized D q family of multifractal dimensions, defined in Balatoni and Rényi (1956) (see also , Rényi (1961)) and later in Hentschel and Procaccia (1983), and independently of them in Grassberger (1983), lends itself also to be considered. D 0 is generally referred to as the capacity dimension, and it coincides with the box-count dimension; hence it is among those mentioned previously and subject of our analysis. Frequent and sudden fluctuation in the spring discharge is needed to produce a rough and irregular hydrograph. The level of roughness measured by fractal dimension corresponds to the intensity-in frequency rather than amplitude-of the fluctuation. The fluctuation results from the irregular arrivals of water from precipitation-and other underground sources when present-and is modified by the fractured karstic aquifer structure. The geometry -characterized by the fractal dimension -of its conduits regulates the modification and shapes the final water discharge. Hence this geometry is responsible for the altered fluctuation, as compared to that of the precipitation. Therefore, the analysis of fractal dimensions of spring discharges may provide a way to obtain information on the hidden subsurface conduits and the flows of water in them in an exploratory and comprehensive manner the physical observations are not capable of. In the present paper, we analyze various fractal dimension estimations and relate them to the complexity of the aquifer and water migration processes in them.

Theory and methods, estimation of fractal dimensions
A geometrically intuitive notion of dimension is an exponent that expresses the scaling of an object's bulk with its size (cf. Theiler 1990). It can also be interpreted as how thickly an object fills part of the space of an integer topological dimension in which it is embedded. Explicitly, for time series, the fractal dimension measures the roughness (or conversely, smoothness) of the paths; the rougher the path, the higher the dimension between 1 and 2. The earliest definition of a non-integer dimension is Hausdorff's; since then, various definitions have become known within this concept, reflected in the estimation methods. Several estimating methods are thematically described in Gneiting and Schlather (2004) and Gneiting et al. (2010). Below we give a short overview of those we use in the present paper.
The Hausdorff Dimension Hausdorff (1918) gives a theoretically rigorous definition of a non-integer dimension of a set. He covers the considered set with balls of a radius at most r , calculates the sum of the d-th power of those radii, and takes its infimum over all coverings of the type. inf coverings r i <r r d i As r → 0, the limsup of this infimum is either infinity or 0, and Hausdorff showed that the transition happens at a unique value dim H as d grows. This transition value dim H is the Hausdorff dimension of the set. For further details, see, e.g., Theiler (1990).

The Box-count Estimator
The infimum taken over all coverings in the Hausdorff dimension makes it difficult, if not impossible, to implement it numerically. Box-counting is the simplest way to overcome this difficulty. Under light regularity conditions, this estimation converges to the Hausdorff dimension. In this method, the path of a time series is initially covered by a single box (square, cube, hypercube, depending on the envelope dimension), which is then divided into smaller boxes of equal size. The method counts how many of them still covers the path. This procedure is to be continued until the box width equals the resolution of the data. For fractals, the growth of the number of boxes is power-like. The standard procedure to obtain the exponent of the power-like growth goes by regression. Regressing the number of boxes on their width on the log-log scale, the regression line's slope, i.e., the regression coefficient, gives the estimated box-count dimension.

Variogram Based Estimators
For a time dependent stochastic process {X t : t ∈ R + } with stationary increments, the V p (t) (semi)variogram function of order p is and its method of moments estimator at lag t = 1 n is (1) When p = 2, we recover the usual variogram, the case p = 1 provides the madogram, while p = 1 2 gives the rodogram. The FD estimator based on the p-th order variogram can be obtained from the regression fit of log(V p )(t) on log t.
The Increment1 estimator is created similarly using p = 1 and in (1) second order differences instead.
In case of p = 1 the relationship between the madogram and fractal dimension is particularly universal and can be explained through the Lipschitz-Hölder heuristics of Mandelbrot (1977), page 304, see also Carvalho and Caetano (2012).
holds for its increments around x. Kahane (1985) Chap. 10 Sec.7 relates Hölder continuity to the Hausdorff dimension. A special case of the result in one dimension, i.e., for real functions over R is that a Hölder continuous function with the supremum of its Hölder exponents β with 0 < β < 1 has a fractal graph, and its Hausdorff fractal dimension is (2 − β). Such a function is irregular, and in general, it cannot be given by a formula. The trajectory of the fractional Brownian motion may serve as an example. Now, let's take the case of the madogram. When regressing log(V 1 )(t) on log t we get the regression coefficient β, this means a power like growth of the increments around t with exponent β, i.e. Hölder continuity with maximal exponent β. Hence the fractal dimension of the observed path is 2 − β. The relationship of FD with the variogram and the rodogram is less straightforward; see Gneiting et al. (2010).
The following assertion is important to us for the interpretation of our results. The sum of two Hölder continuous functions with (maximal) exponents p and q is again Hölder continuous with (maximal) exponents min( p, q). Hence, the sum of two independent fractal processes' fractal dimension is equal to the larger fractal dimension.
The mentioned estimators are available in the [R package 'fractaldim'], 1 whereas the generalized fractal dimension estimators are available in the [R package 'fractal']. 2 As is well-known (Adler and Taylor 2007), the Brownian paths are nowhere differentiable, not even Lipschitz continuous. However, they are Hölder continuous in all t with any of the exponents 0 < β < 1 2 and their Hausdorff dimension is 3 2 . The paths of the fractional Brownian motion of Hurst exponent H are Hölder continuous in all t with any 0 < β < H while their Hausdorff dimension is 2 − H , showing, that the bound in Kahane's theorem is sharp.
The fractal dimension of independent identically distributed (i.i.d) series is a specific issue, since i.i.d random variables indexed by continuous time constitute a nonmeasurable process; hence they are not a usual subject for probabilistic analysis. Nevertheless, for i.i.d. random variables with a continuous distribution, the "limiting graph" as the index set tends to continuum, covers a full segment of the plane and therefore its fractal dimension is 2.

The Rényi or generalized D q Family of Multifractal Dimensions
In computing the box-counting dimension of a set, no provision is made for weighting the box count according to how many points are inside. The generalized (multifractal) dimension D q does take into account the number of points in the box. For a monofractal or uniform fractal, D q does not vary with q. For a non-uniform fractal, i.e., multifractal, the variation of D q with q quantifies the non-uniformity. Hence the set of D q -s is also regarded as multifractal dimensions. For q = 2 D q is called the correlation dimension, while for q = 1 the information dimension. Let us note that q = 0 corresponds to the plain box-counting dimension defined above and is often also referred to as the capacity dimension.
In order to analyze the relationship between the multifractal dimensions and the fractal dimensions previously mentioned, the concept of phase space embedding will require some elaboration (for more details, see, e.g., Theiler (1990), Addison (1997), Kantz and Schreiber (2004), and references therein.) Instead of studying the temporal course of spring discharge through the hydrograph, it is possible to regard the phenomenon as the output of a dynamical system driven by a random noise generator and then consider its time evolution in some phase space. Given the system's mentioned non-linearity, irregularity naturally presents itself in the output, often creating a multifractal object or attractor in the phase space. Although the state space description is highly beneficial from the theoretical perspective, an immediate question to face in practice is the state space's unknown dimension. Several embedding theorems clarify the proper phase space dimension choice when the geometrical object formed by the phase vectors in the phase space is equivalent to the trajectory, the hydrograph. Under quite general circumstances, this dimension is larger than twice the box-counting dimension (Addison 1997). When the dimension of the underlying phase space is not known (Addison 1997), advocates a practical method to calculate information and correlation dimensions, the focus of our interest, in embedding spaces of successively larger dimensions. The estimated dimension (the slope in the log-log plot) will initially increase with the growing embedding dimension, reaching a limiting value when the embedding space is large enough. When the increase of the embedding dimension stops growing the slope, the obtained value can be regarded as the estimation of the true correlation dimension. The paper follows this methodology.
There are other methods readily available for fractal dimension estimation. However, those have not proven to be useful for the present study. E.g., frequency domain and wavelet-based estimators unanimously gave 2 for the fractal dimension for quite a few springs. That is obviously anomalous behavior.

Geological and hydrological settings
Aggtelek Karst is located on the Hungarian Slovakian borderland, in Gömör-Torna Karst, which belongs to the Inner West Carpathians. Gömör-Torna Karst is divided into Slovak Karst (northern part) and Aggtelek Karst (southern part) (see Fig. 1). The area is a geographically homogeneous region extending over 600 km 2 . The Hungarian side consists of East-West orientated limestone plateaus, having a surface area of about 200 km 2 and a height of 400-600 meters a.s.l. Deep valleys are incised between the plateaus, with their lowest points being at an elevation of 150-260 meters. On the other side of the state border, the highest plateaus of the Slovak Karst have a height up to 800 meters, and to the north, the mountain range is even higher (Haas 2012). More than 1000 caves of this karst area were listed as part of UNESCO's World Cultural and Natural Heritage in 1995.
The diverse geology of the Aggtelek Karst inspired numerous studies, and as a result, it became Hungary's best-explored karst area (see, e.g., Láng (1955), Jakucs (1977), Hevesi (1991), Telbisz et al. (2020), and references therein). The main mass of the mountain is built up from Triassic carbonate rocks, which were deposited on the carbonate platforms of the Neotethys Ocean (see Fig. 2). The Lower Triassic layers are mainly aquicludes; they contain clay and sandstones, covering about 62 km 2 . Above them, the Middle and Upper Triassic layers consist of well karstifiable limestones (Wetterstein Formation is the widest formation), and to a lesser extent, of less karstifiable dolomites with an area of about 105 km 2 . Sand and gravel layers from the Miocene age are much smaller in the southern part of the area. In the Miocene, this area was flooded by the Pannonian Sea. The Sea later gradually decreased, became a lake, and lacustrine sediments were deposited on the carbonated area (Less 2000). In the second part of this period, karstification began under subtropical climate conditions. In the Pliocene, the north part of the Karst area started to uplift, and due to the tectonic forces, the present-day dominant north-south sloping surface was formed. Pliocene sediments (clay, sand, gravel) cover an area of approximately 35 km 2 . The Pleistocene  level is represented by the 0.2-1.0 m thick clayey soil on the plateaus' top. In the valleys, the Holocene alluvium consists of gravel and sand. According to Grill et al. (1984), the mountain structure can be described by assuming the existence of overthrust folds. Four different overfields can be distinguished. The plateaus are situated above pseudo-synclines, where they remained untouched between the uplifted tectonic zones, which brought to the surface the Lower-Triassic layers.
The mean cumulated annual precipitation is around 600 mm. The greater part of the precipitation infiltrates the dolina-covered plateaus and arrives at the springs via the fracture-system network. The biggest karst springs rise along the foot of the hills or at the end of closed valleys and the contact of karst and non-karst areas, while smaller ones can be found at higher altitudes. The area has a humid continental climate, with a long summer and a substantial effect on the mountain's climate due to the proximity of the Carpathians. The mean annual air temperature is 9.1 • C. In a karstic environment, the air temperature in caves and conduits is almost constant throughout the year and is nearly equal to the mean annual surface air temperature. The area is characterized by 120-130 days of freezing temperature and 40 snow-covered days in a year. The average annual cumulated snow depth is 0.5 m.

Data description
Measurement of the daily spring discharge of the most significant karstic springs in the Aggtelek Karst was carried out by the Hungarian Water Resources Research Plc. from the early 1960s to the 1990s, see Maucha (1998). We could select 12 springs with non-missing daily observations for 6205 days in the period of 1 March 1974 to 23 February 1993. The basic descriptive statistics of these decades of uninterrupted hydrological data series are shown in Table 1. The surface area of the catchments is spread over various scales ranging from 0.5 km 2 to 34.4 km 2 . The springs can be found at altitudes between 165-268 m a.s.l. The primary source of water in the springs is the local karstic infiltration from the surface. The average and maximum water temperature of the springs is around 10 • C, just as the mean annual air temperature in Hungary mentioned in the previous chapter. However, as Table 1 shows, significantly higher average and maximum temperatures were observed at several springs. The most extreme is Melegvíz ("Warm Water", coded as Mel) spring with an average water temperature of 17.8 • C, and maximal water temperature of 19.4 • C. The springs are located within a relatively small area, close to each other. Nevertheless, significant differences of scale can be observed in their discharges: the multiannual averages vary between 255 and 14509 m 3 /day, see Table 1. The lowest average discharge was measured at Mel spring, while the two springs having the largest catchment areas, Jósva (Jos) and Nagy-Tohonya (NaT), have the highest averages. The appearance of extreme discharges can be characterized by comparing the usual discharge -Q med represented by the median discharge value in order to circumvent the biasing effect of high values in the mean -with the Q max maximal ones by the Q max /Q med ratio. The ratio varies between 3.1 and 160 and is the smallest for springs Mel, Tap, NaT, Jos (being 3.1, 15.4, 19.2, and 28.9, respectively) having higher than usual maximum temperature. Again, this may be the effect of the flow system mentioned earlier, which provides a steady water supply in these springs and does not allow for low values in the denominator, keeping the ratio low.
Cave systems of varying lengths are to be found in connection with some of the springs (Kordos 1984;Székely 2003). At Jos and Mel springs, water comes to the surface from the 25 km long Baradla-Domica cave system, the second-longest cave in Hungary. Kom spring is connected to Béke cave, of length approximately 7.2 km. Szabadság cave of approx. 3.2 km belongs to the aquifer of Kec spring. NaT spring transports water from the 1.4 km long Kossuth cave to the surface. The Vass Imre cave of length approx. 1.1 km can be found at KiT spring. Caves of much shorter confirmed lengths -typically a few hundred meters, as cave researchers put it by current knowledge -lie behind Cso (Csörgő cave), Kop (Kopolya cave), Bol (Róka cave) springs.

Results of fractal dimension estimations
Fractal dimensions of spring discharges are computed as described in Sect. 2. Their values are given in the Appendix, in Table 2.
The various methods emphasize different aspects of fractality or scaling, and as a result, the obtained values differ significantly (Fig. 3). It is not possible to single out one particular method that would sufficiently comprise the information provided by all of them.
Concerning the Rényi D q family, the embedding dimensions turn out to vary widely for the springs. The obtained information dimensions vary between 1.74 and 3.57, and the correlation dimensions between 1.83 and 3.52. The magnitude of values is not surprising since these are FDs of objects in the phase space that are only equivalent in the sense of Sect. 2 to that of the trajectories.
So, mixing them with the FDs of the trajectories in the analysis would not be meaningful. However, in order to have a visual impression, we display the box-count (D 0 or capacity in the multifractal terminology) and the linearly rescaled -for the sake of better visual comparability -D 1 or information dimensions on the lower plot of Fig. 3. It is apparent that the values of the estimated D 1 information dimension have When computing correlation dimension D 2 for the springs, the procedure did not converge for 2 of them, Vec and Kas, meaning that the values kept steadily growing with the increase of the embedding dimension up to 16. For three other springs: NaT, Mel, and Tap, the FD stabilized over embedding dimensions 9, 10, and 12, which seems to be suspiciously high. On the other hand, for Jos, Kec, KiT, and Kom, the correlation dimension almost coincides with the information dimension-just as is typical for monofractals. Since at least five springs do not have reliable correlation dimension values, it makes more sense to go on with the analysis and interpretation of the results without the values of D 2 .
Grouping of springs may help in explaining the variability of FD values across springs. As the box-count dimension of spring hydrographs shows a different pattern across springs from the rest, two different groupings were made and represented by their dendrograms in Fig. 5. The upper dendrogram is based on box-count FD values. For this single variable, clustering reduces to a ranking by the magnitude of values. This practice may not be regarded as a proper cluster analysis-even though the mathematical procedure is meaningful in this case, too. It is a standard procedure in the literature (mainly medical or biological statistics) for obtaining the segmentation of non-homogeneous observations. The other clustering, represented by the lower dendrogram in Fig. 5, has been made by all p-variogram-based FDs. The calculations were done by Ward's method, and since the p-variogram related FDs are correlated, Mahalanobis' distance has been used. To help understand and interpret the results of cluster analysis the hydrographs of the groups are displayed in Fig. 6 in accordance with the upper dendrogram in Fig. 5, i.e., the box-count dimension.
It is well observable in Fig. 6 that roughness, i.e., the intensity of the hydrograph's fluctuation, is similar within groups. The hydrograph is visibly the smoothest for the standalone Jos spring and the roughest for the second group. This is well in line with the fact that the box-count dimension estimates the Hausdorff dimension, which is just the roughness's theoretical measure. The difference between the two subgroups of the first group does not seem perceptible from the hydrographs by naked eyes.

Similar approaches in the literature
River networks have long been noted for their fractal and multifractal nature (La Barbera and Rosso 1989). In this context, generalized fractal dimensions have been assessed through different methods. It was shown that the multifractal approach char-acterizes the river network and also the flood hydrograph as the basin response to heavy rainfall events (Gaudio et al. 2004). Similarly, the karst media (matrix rock porosity, fracture network, and karst conduit network) are also fractal . A recent review of the topic with rich references to the literature is Pardo-Igúzquiza et al. (2019). However, the karst media and spring hydrographs are not directly comparable to the river networks and hydrographs, although similar methodologies in the analysis may be applied. The functioning of the karst system is a widely studied field (cf., e.g., Bauer and Tóth (2015)). Veneziano and Essiam (2003) have shown that for an aquifer with a multifractal hydraulic conductivity field, the resulting hydraulic gradient and specific discharge are also multifractal fields; see also (Florea 2001;Skoglund and Lauritzen 2011). Hence, an analysis focusing on the fractal or multifractal characteristics of the discharge of springs may shed light on certain properties of the karst media that would otherwise be difficult to infer directly. In some cases, eventual inconsistencies with existing geological knowledge may point to insufficient information in certain aspects and thus inspire and direct further researches, e.g., further exploration of caves, etc. In Robledo-Ardila et al. (2014), four springs are taken in different carbonate massifs, and their fractal property is analyzed in the above-mentioned spirit. Their results show that each spring's temporal distribution of flow has a low and varying fractal dimension. When translated to the graph of the discharge hydrograph as a planar object, it turns out to be in the range of 1.3-1.45, just as our springs with a simple flow system, i.e., without mixing with a hydrothermal or any other component. They claim that, on the whole, the fractal dimension is low. However, there are temporal periods when the fractal dimension is high, indicating multifractal properties. We do not consider temporal periods but establish the multifractal property for some of the considered springs by the moment scaling. Hergarten and Birk (2007) present an analytic theory for the recession of spring hydrographs after precipitation events, focusing on short time scales. In order to validate their theoretical considerations, they analyze two recession curves of 10 days length of the spring Gallusquelle. They find 1.4 as the box-count fractal dimension, which is similar to our findings; however, the two approaches are entirely different. Our 6502 days long observations contain hundreds of recession curves, allowing for a recession curve analysis only of exemplary nature. Therefore, this approach, followed also in Fiorillo (2014), is not followed here.

Mono-or multifractal
Since we intend to analyze fractal dimension estimations, a detailed multifractal analysis is very much beyond the scope of the present paper. Nevertheless, we present the multifractal spectra of the springs in the Appendix, in Fig. 9, computed by the [R package 'MFDFA'] 3 (its code is identical with the corresponding package of MATLAB). Those can be compared to the fractional Brownian motions' multifractal spectra, also given in the Appendix, in Fig. 10. It is noteworthy that the peak of f (α) occurs at or very close to 1 + H = 3 − D 0 . For this reason, we draw a vertical line at three minus the box-count dimension: 3 − D 0 when displaying the multifractal spectrum of the spring hydrographs in Fig. 10 in the Appendix. The picture does not change much with the changing of the scale and the power orders (those are input parameters of the computational package). For the two springs with no convergent correlation dimension, a clear peak happen to be missing in the graph, and the vertical line does not fit for the location of the peak of a hypothetical curve fitting to the points. This discrepancy may give a hint of why the correlation dimension estimations do not converge for these two springs. Beyond that, however, no apparent pattern seems to explain the very high embedding dimensions. Again, we do not go into the analysis of this question.
Moment scaling, presented in the Appendix in Fig. 12, indicates that most hydrographs are monofractals rather than genuine multifractals (cf. Schertzer and Lovejoy (2011), section 6.3). Hence, a single number for the fractal dimension is informative. That may also explain why box count and information dimensions' relation is linear. Theoretically, for monofractals, they should coincide. In our case, they differ by roughly 1, precisely the difference in the embedding dimension (remember the graphs are shifted in Fig. 4 for better visibility only). Therefore, they similarly characterize the springs, and as a result, information FD does not provide valuable new information but strengthen that of the box-count (capacity) dimension, which is also important.

Uncertainty of the FD estimation
Addressing the computed values' uncertainty is essential, but there is currently no available distribution theory for fractal dimension estimators. The few existing uncertainty estimates are model specific, and those models have to be fitted beforehand. In the present study, we do not fit any model to the hydrographs but address the issue by a general approach. We only try to obtain an idea about uncertainty when the true fractal dimension is known. We consider a simulation from fractional Brownian motions of the Hurst coefficient H and compute the FD estimates' deviation from the theoretical value 2 − H . Those estimated FD values are displayed in Fig. 13 in the Appendix, and the average values and standard deviations are presented in the Appendix in Table 3. It could be claimed with high certainty that the differences between the FDs of springs are well above the random fluctuations caused by error or noise in the observations, were those values valid for the hydrographs.

FD of precipitation data
The discharges of karstic springs are influenced primarily by precipitation. So, even though we do not intend to model the discharge-precipitation relationship, it is worth determining the FD of precipitation data for a starting point. Because of the relatively small study area, the discharges of springs are influenced by very similar precipitation in terms of both quantity and character (rain, snow, etc.). Thus, one monitoring station's daily registered precipitation data (Jósvafő) located in the region is sufficiently representative for our purposes. Besides, for the given time period, no other uninterrupted monitoring data for precipitation is available. The average fractal dimension computed for these data by the paper's methods is 1.81, a high value. The standard deviation of 0.08, so the values are close to each other. In Gneiting et al. (2010), the authors argue for using the madogram as the most reliable fractal dimension estimation. In the given case, the madogram value is 1.83, underlining again that the graph of precipitation has high FD and agrees well with FDs of i.i.d. series. The non-zero values of daily precipitation are known to be close to probabilistic independence, and the zeros have no dominating role in shaping the graph's smoothness. This circumstance manifests in the fact that the graph fills almost completely a segment of the plane, i.e., the mentioned fractal dimension is close to 2.

Effect of infiltration
The process of infiltration and migration of water in the karstic environment holds back a significant amount of water arrived in the precipitation and, as a result, creates much smoother hydrographs for springs than precipitation data. Consequently, the hydrographs have much smaller FD than precipitation. However, the smoothing effect is not uniform in the study area. It depends on the structure of the spring's aquifer, particularly the conduits' geometry, and on the mixing of water from other sources or subsurface flows if they exist. Therefore, despite the similarities in the territorial distribution of precipitation and the area's large-scale geological structure, the springs' FDs vary significantly. E.g., madogram FD values vary from 1.13 to 1.61, and boxcount values vary from 1.31 to 1.69, representing the differences between the springs.
Effect of additional sources of water Examining the graphs of FDs in Fig. 3, a striking difference between the boxcount dimension and the variogram-based dimensions catches the eye. Therefore, the two types of fractal dimensions are likely to represent different characters of the aquifers of springs. The results of cluster analysis may help to explore this difference. Consider first the box-count dimension and the clustering of springs according to it as displayed in the upper graph of Fig. 5. The three rightmost springs: Mel, NaT, Tap in the dendrogram have a significantly higher maximum water temperature than the others (see Table 1). The elevated temperature can be explained by the fact that another water source controls these springs' outflow beyond the immediate effect of infiltration. Probably, intermediate flow systems bring the water of somewhat higher temperature and combine with the water arriving from infiltration (Goldscheider et al. 2010). The presence of additional water source is supported, e.g., by tritium examinations of the waters of Nagy-Tohonya spring, where 1 to 4 years old water is found to reach the surface (Deák and Dénes 1981), postulating a complex system with more components. The flow system is more or less independent of infiltration; therefore, its effect can particularly well be observed at low discharge when supply from infiltration is low and outflow temperature is higher. Furthermore, their maximum discharges are low compared to the usual discharges; the Q max /Q med ratios are 3.0, 19.2, and 15.4, respectively. That suggests that another, more steady, subsurface flow system mixes with the infiltrating water, as described in Sect. 4. The superposed effect of two water sources, -the infiltration and the intermediate flow-creates a complex fluctuation pattern in the hydrograph and results in higher box-count FD. Indeed the box-count FD is 1.69, 1.54, 1.49 for these springs. The bulk of the springs is arranged in two subgroups of the largest group at a relatively low level of separation. Remarkably, one of the subgroups (Kas, Bol, Kit, Kop) has the smallest difference between the maximum and minimum temperatures of their water, 1.4, 1.3, 1.3, 1.2 • C, respectively. The temperature differences for the other subgroup (Kom, Cso, Vec, Kec) are 2.8, 3.2, 1.8, 2.2, • C, respectively. We do not know any particular geological or hydrogeological character that would justify just this separation. However, in the first mentioned subgroup, the maximum water temperature adapts slightly more to the region's general cave temperature than the second subgroup (except for Kec). That suggests, the water either spends more time in the conduits of the aquifers or the conduits' structure is finer to accelerate heat exchange.
The leftmost (in the dendrogram) standalone spring Jos is in a unique situation. It has the largest aquifer and collects the highest amount of water that manifests in very high discharges. Consequently, its hydrograph has very high peaks related mostly to large precipitation events (Fig. 6). Compared to those peaks, the regular regime's fluctuation is small, and therefore, the hydrograph's geometry is smoother at this resolution. The high peaks and small fluctuation create a uniquely low box-count FD of 1.31 for this spring. At the same time, both the higher maximum temperature and the lower Q max /Q med ratio (28.9) would associate this spring to the second, rightmost group. Relations of the box-count and madogram FDs with the water temperatures and the discharge ratios are displayed in the scatterplots given in the Appendix, Fig.  8. The hydrogeological situation is that its aquifer overlaps with that of Mel spring; therefore, the intermediate flow influences it, but the huge amount of water supply from infiltration suppresses its effect.
Effect of karstification It lends itself to study the dependence of the FDs on karstification, as introduced in Sect. 3. When precipitation reaches the area, under suitable conditions, the cave system gives way to a quick outflow of a large amount of water, creating very high peaks in the hydrograph. On the other hand, the karst's fractured system holds back part of the water and discharges it with a non-uniform delay, dependent on the structure. These processes are similar, and the created water discharges decay with different exponents over time (Maucha 2005); their superposition forms the spring discharge. This superposition creates variable similarities when varying scales are considered, typical of how fractal and multifractal processes are constructed by nature. As noted in Sect. 2, the Lipschitz-Hölder heuristics shifts the fractal dimension of superposition of functions towards the highest FD of the components. Therefore, more variable delayed discharge of water creates more complex fluctuation patterns in the hydrographs reflected by increased FDs.
When the springs are clustered according to the V p family-related FDs: variogram, madogram, rodogram, and Increment1 based ones, three groups are clearly distinguished, as shown in the lower dendrogram in Fig. 5. The springs of the first group: Mel, Jos, NaT, Kec, Kom, and Cso all have large cave systems in their aquifer, an indicator that the level of karstification of their aquifer is higher. Note here that the known length of the cave in the aquifer Cso is smaller than that of the others. The springs in this group have high FDs, e.g., madogram FDs vary from 1.29 to 1.61. All springs of the other group have aquifers with a low karstification level. Their discharges present less complex fluctuation patterns, and the result is a smoother hydrograph reflected in lower fractal dimension changing between 1.13 and 1.21.
All four types of FDs contain valuable information because leaving out any of them results in a similar but not identical grouping. When omitting any of the springs from the analysis, the grouping of the remainder does not change. Hence the grouping is stable and reflects the known geological structure influencing the springs.

Combined effects
Ideally, the two types of fractal dimension estimations considered in this section are meant to estimate the same quantity. In the limit of infinitely frequent observationspractical nonsense-they are supposed to give the same values. The fact that the different procedures still provide different values may either be due to the data's finite resolution or the violation of the ideal conditions. As a result, the different methodologies emphasize different effects. However, it does not mean they would be independent characteristics, and they would ignore the other effect. The effects' combined appearance is reflected in the lower dendrogram of Fig. 5. The group of six springs with the highest karstification is divided into two subgroups. That division is nothing to do with the karstification level. Instead, the subgroup on the right-hand side consists of NaT and Mel, where the flow mixing effect is strongly present, and (Jos) where it exists but weak, while the mixing effect is lacking in the left subgroup (Kom, Kec, and Cso).
Calculation of fractal dimensions was essential for obtaining the grouping by the mentioned principles. Clustering by the hydrographs themselves would lead to the result presented in Fig. 7. Remarkable that Jos and NaT separated at a very high level from the rest happen to create flash floods while the others do not do. Beyond this, however, it does not seem easy to interpret this result geologically.

Conclusions
Discharge data registered daily over more than two decades for 12 springs of the Aggtelek karst were analyzed with a focus on the fractal properties of the graphs of the registered time series, i.e., the hydrographs. The springs are located in a geographically compact and thoroughly explored area of a similar lithological structure.
The infiltration and migration of water in the fractured karst make the temporal course of the amount of discharged water in springs much smoother than the course of precipitation reflected in a hydrograph of much lower FD. The estimations of FD capture different aspects of fractal nature. The p-variogram based estimations may be associated with the karstification level. At the same time, box-count and information dimensions may reflect the mixing of discharge sources like regional flow components acting independently of infiltration. The analysis of FDs can therefore strengthen or broaden our knowledge on the aquifers of springs. An FD not in line with our prior knowledge may indicate that the aquifer is more complex, including not well-explored areas or multiple, not well-explored sources that also provide water to the springs. These are critical aspects of aquifer vulnerability as contamination in an area not known to be associated with the aquifer may pollute its water. Similarly, contamination may arrive with subsurface flows. The lack of their adequate exploration precludes the implementation of suitable protection measures and the aquifer's water resources may be spoiled (Fehér et al. 2016). Having this in mind, the FDs' presented analysis has a potential to be an invaluable environmental research tool.
Further work may concern the proper modeling of the discharges. As theory says, linear models such as diffusions described by stochastic differential equations invariably create processes with paths of FD value 1.5. So, modeling has to go beyond those and to find the proper class of models is undoubtedly a challenge.