Influence of initial heterogeneities and recharge limitations on the evolution of aperture distributions in carbonate aquifers

Karst aquifers evolve where the dissolution of soluble rocks causes the enlargement of discrete pathways along fractures or bedding planes, thus creating highly conductive solution conduits. To identify general interrelations between hydrogeological conditions and the properties of the evolving conduit systems the aperture-size frequency distributions resulting from generic models of conduit evolution are analysed. For this purpose, a processbased numerical model coupling flow and rock dissolution is employed. Initial protoconduits are represented by tubes with log-normally distributed aperture sizes with a mean μ0 = 0.5 mm for the logarithm of the diameters. Apertures are spatially uncorrelated and widen up to the metre range due to dissolution by chemically aggressive waters. Several examples of conduit development are examined focussing on influences of the initial heterogeneity and the available amount of recharge. If the available recharge is sufficiently high the evolving conduits compete for flow and those with large apertures and high hydraulic gradients attract more and more water. As a consequence, the positive feedback between increasing flow and dissolution causes the breakthrough of a conduit pathway connecting the recharge and discharge sides of the modelling domain. Under these competitive flow conditions dynamically stable bimodal aperture distributions are found to evolve, i.e. a certain percentage of tubes continues to be enlarged while the remaining tubes stay small-sized. The percentage of strongly widened tubes is found to be independent of the breakthrough time and decreases with increasing heterogeneity of the initial apertures and decreasing amount of available water. If the competition for flow is suppressed because the availability of water is strongly limited breakthrough of a conduit pathway is inhibited and Correspondence to: B. Hubinger (bernhard.hubinger@uni-graz.at) the conduit pathways widen very slowly. The resulting aperture distributions are found to be unimodal covering some orders of magnitudes in size. Under these suppressed flow conditions the entire range of apertures continues to be enlarged. Hence, the number of tubes reaching aperture sizes in the order of centimetres or decimetres continues to increase with time and in the long term may exceed the number of large-sized tubes evolving under competitive flow conditions. This suggests that conduit development under suppressed flow conditions may significantly enhance the permeability of the formation, e.g. in deep-seated carbonate settings.


Introduction
Karst aquifers are important drinking water resources but highly vulnerable to contamination due to rapid transport of pollutants through solution conduits.Solution conduits evolve as a consequence of dissolution processes enlarging discrete pathways along fractures or bedding planes.With increasing aperture size the flow resistance decreases, thus accomplishing higher flow rates.In turn, higher flow rates enhance dissolutional enlargement.Thus, there is a positive feedback between dissolution and flow, which finally leads to a breakthrough event, i.e. strong increases of dissolution rates and discharge within a short time period.
A quantitative understanding of this feedback mechanism was provided by numerical models (e.g.Dreybrodt, 1990;Palmer, 1991) simulating the evolution of mature solution conduits from intially narrow protoconduits.To support the assessment of vulnerability of karst aquifers these models may be employed to identify parameters that control the evolution of conduit systems and thus the patterns of fast flow pathways.
A number of modelling studies (e.g.Howard and Groves, 1995;Siemers and Dreybrodt, 1998;Dreybrodt and Siemers, Published by Copernicus Publications on behalf of the European Geosciences Union. B. Hubinger and S. Birk: Influence of hydrogeological conditions on aperture distributions 2000; Dreybrodt et al., 2005;Kaufmann and Braun, 1999) suggest that evolving conduit systems range between two extremes, namely simple systems with only few, but strongly enlarged conduits as opposed to more complex systems with many interconnected, but less strongly dilated pathways.To describe these two different end-members Worthington (2009) proposed the terms "macrokarst" and "microkarst".
A more quantitative approach to distinguish different types of conduit systems was introduced by Bloomfield et al. (2005).
These authors plotted the cumulative aperture distribution resulting from numerical simulations of conduit evolution.They generally found bimodal aperture distributions, where a certain number of large conduits evolve while all other apertures stay close to their initial values.They only consider Darcian flow and assume aperture growth rates being proportional to volumetric flow rates to some given power in the range of 0.1 to 0.8.The number of large-aperture conduits was found to be controlled by the chosen exponent of the generic rate law used for describing dissolutional widening and also by intial aperture heterogeneity.Rehrl et al. (2008Rehrl et al. ( , 2010) ) used cumulative aperture distributions to analyse conduit evolution in artesian gypsum karst settings of the Western Ukraine.Bimodality was found to depend on the chemical saturation of the inflow but the effect of the hydraulic conductivity of the rock formation incorporating the conduit network was even more pronounced.With low hydraulic conductivity a smooth transition from small to large apertures rather than a bimodal aperture distribution was found to emerge.It was thus suggested that the limitation of flow imposed by low hydraulic conductivity provides an important control on the geometric properties of the evolving conduit systems.
The main objective of this work is to provide a closer examination of the potential effects of limited flow rates on aperture distributions evolving in carbonate rocks.Further, the effect of the initial heterogeneity on aperture distributions is considered to provide a link between results by Bloomfield et al. (2005), who used a generic dissolution rate law, and earlier studies (e.g.Howard and Groves, 1995;Siemers and Dreybrodt, 1998) that accounted for carbonate dissolution kinetics but did not evaluate aperture distributions.Thus, the numerical simulations aim to address effects of both geological (initial heterogeneity) and hydrological (limitation of flow) conditions on conduit evolution in carbonate rocks.

Modelling approach and setup
The tube network flow model implemented in the modelling tool CAVE (Carbonate Aquifer Void Evolution; e.g.Liedl et al., 2003;Birk et al., 2006) is used to represent networks of protoconduits that are potentially enlarged due to rock dissolution.According to Kirchhoff's law (Clemens, 1998;Horlacher and Lüdecke, 1992) total inflow and total outflow have to balance at each node of the tube network.Depending on the average (skalar) flow velocity v (ms −1 ) within a tube of diameter a (m), the dimensionless Reynolds number N Re given by (cf.Dreybrodt, 1988;White, 1988) is calculated, where ρ w (kgm −3 ) is the density of water and η w (kgm −1 s −1 ) the dynamic viscosity of water.Both ρ w and η w are calculated according to Weast (1979) for a temperature of 10 • C. N Re is used to distinguish between laminar and turbulent flow.Flow switches from laminar to turbulent flow for N Re ≥2000 and from turbulent to laminar flow for N Re ≤ 1000.Discharge through single tubes is calculated using the Hagen-Poiseuille-equation (White, 1988) for laminar flow conditions.In case of turbulent flow, the Darcy-Weisbach-equation (White, 1988) and Colebrook-White-equation (Dreybrodt, 1988) are employed (Horlacher and Lüdecke, 1992;Hückinghaus, 1998).
Once flow rates in the tube network are known, rock dissolution is calculated.In general, the dissolution rate F (molm −2 s −1 ) can be described by where c (molm −3 ) is the calcium concentration in the water, c eq (molm −3 ) is the calcium equilibrium concentration with respect to limestone, k (molm −2 s −1 ) is a reaction rate constant, and n is a dimensionless exponent.If diffusion is the rate limiting process k = N Sh Dc eq a (Dreybrodt, 1988;Bauer, 2002;cf. Bauer et al., 2003;Dreybrodt et al., 2005), where D (m 2 s −1 ) is the diffusion coefficient of water and N Sh is the dimensionless Sheerwood Number.For laminar flow we assume N Sh = 3.66 (Beek and Muttzall, 1975) and for turbulent flow N Sh is calculated according to Gnielinski (1976).If c is far from c eq n equals unity (Buhmann and Dreybrodt, 1985a,b;Dreybrodt, 1988;Dreybrodt et al., 1996).Experimental findings (Svensson and Dreybrodt, 1992;Eisenlohr et al., 1999) suggest that dissolution rates drop dramatically if c exceeds a certain switch concentration c sw (molm −3 ), which is close to c eq .An overview of calcite dissolution kinetics is given by e.g.Kaufmann and Dreybrodt (2007) and Dreybrodt et al. (2005).
The calcium concentration is calculated at each node assuming complete mixing of inflowing waters.
The resulting concentration is assigned to water entering downgradient tubes.
Each tube is subdivided into 25 equidistant segments.Depending on dissolution rate and input concentration, approximations given by Bauer (2002) are used to obtain the calcium concentration at the outlet of each segment, which is passed to the next segment in flow direction.According to the law of mass conservation, the mass dissolved along a small distance equals the mass difference of calcium flowing into and out of the segment.Thus, further considering the rock area exposed to water and calcium concentration c Ca (molar density; molm −3 ) in the rock, the carbonate mass dissolved from the wall is transformed to a change of aperture size.Finally, a new hydraulically equivalent aperture size is calculated for each tube.
Flow rates and apertures are updated every 10 a and for the long-term analyses described in Sect.4.2.2 every 1 ka.In preceding tests, this time discretisations were found to represent the changes with sufficient accuracy for the simulations presented within this work.Table 1 lists parameter values used for the calculation of dissolutional enlargement.They refer to limestone and are geared to Dreybrodt et al. (2005).

Model scenarios
The model scenarios consider a two dimensional area of 200 m length and width, respectively.The modelling domain represents a bedding plane section in a fractured limestone aquifer (Fig. 1), where initially narrow discontinuities (protoconduits) are recharged with chemically aggressive water and therefore can be enlarged by the dissolution processes.The protoconduit network is represented by 722 10 m long water-filled tubes, which are hydraulically interconnected at 400 nodes.Water enters the system at nodes located at the left hand side.Discharge is restricted to nodes at the opposite side, whereas the remaining sides are no flow boundaries.Thus, the model setup corresponds to the generic setting considered by Bloomfield et al. (2005).
Initial protoconduit apertures are spatially uncorrelated and log-normally distributed (cf.e.g.Siemers and Dreybrodt, 1998; Romanov et al., 2002Romanov et al., , 2003;;Gabrovšek et al., 2004).The mean of the logarithm of the apertures is set to µ 0 = 0.5 mm, and the standard deviation σ 0 is varied to account for different initial heterogeneity.For models representing slightly heterogeneous networks σ 0 = 0.1 mm and for more heterogeneous ones σ 0 = 0.5 mm is used.This results in an initial degree of heterogeneity of σ 0 µ 0 = 0.2 and 1, respectively.For each model 30 realisations are performed.The random numbers defining the initial apertures are taken from the same seed for the corresponding realisations of the different scenarios.Thus, the initial location of smaller and larger apertures of a certain realisation is independent from σ 0 but differs for each of the 30 realisations.To ease appraisal of results the same model realisations are reused for each of the scenarios with varied maximum recharge rates described in the following.At early evolutionary stages constant-head conditions with a total hydraulic head difference of h = 5 m are applied at the recharge and discharge sides of the networks.To account for limited availability of recharge inherent in any type of hydrogeological environment, a constant-flow boundary is employed at nodes located at the recharge side at later stages.The total inflow rate Q in (m 3 s −1 ) into the system is limited by a predefined maximum which is varied for different scenarios.The maximum inflow rate Q max (m 3 s −1 ) is split evenly to all i nodes located at the recharge side from which each can receive a maximum inflow rate of q max (m 3 s −1 ), whereby Q in ≤ Q max = i •q max .The constant-head boundary condition at a node is replaced by a constant flow rate boundary condition as soon as the inflow rate reaches q max .Backswitches of boundary conditions from limited inflow rates to constant heads do not occur in the simulations.

Results
If not otherwise stated, the analyses of the model results refer to aperture-size frequency distributions using 10 logarithmic equidistant bins per order of magnitude in the aperture range.From now we will use the term aperture distributions for simplification.References to animations (Anim.)correspond to visualisations given in the Supplement.

Influence of protoconduit heterogeneities
The following subsections examine influences of initial protoconduit heterogeneities on conduit systems evolving in carbonate rocks.

Conduit patterns
A typical example of a conduit pattern at a mature stage of karstification is shown in Fig. 2. The illustration shows one realisation with high initial protoconduit heterogeneity and intermediate total maximum inflow rate.The corresponding temporal evolution of the conduit pattern between 14 ka (shortly before breakthrough) and 30 ka (mature stage) is shown in Animation 1.
In this work the terms "breakthrough" and "breakthrough time" refer to the situation when first-order dissolution kinetics becomes active along an entire pathway.This situation arises as a result of the positive feedback between increasing flow and dissolution rates.Because of the small apertures flow rates are initially low and the calcium concentration approaches the switch concentration still within the tubes connected to the recharge side.Therefore, the downgradient tubes are slowly enlarged following a higher-order dissolution rate law.As a consequence, the discharge through the tube network increases and the region of first-order kinetics penetrates downgradient.When the more efficient first-order kinetics reaches the discharge side the outlet of the tube network is rapidly widenend, thus causing a strong increase of flow and dissolution rates.This feedback mechanism is frequently described in modelling studies on karst evolution (e.g.Dreybrodt et al., 2005) and is also operating in the model scenario considered here.Yet the feedback is finally switched off because the discharge cannot exceed the pre-defined available maximum recharge.
Generally, a main branch is found to develop which usually fans out at the discharge side (Fig. 2).Typically several additional conduits start to evolve at the recharge side but join the main branch soon.This is explained by the re-arrangement of the hydraulic head distribution after the breakthrough (see e.g.Dreybrodt et al., 2005): the breakthrough and the subsequent flow limitation at the corresponding input node causes a head drop along the highly conductive pathway and thus steep hydraulic gradients between the other input nodes and that pathway evolve.As a consequence, there is a rapid development of tributaries that connect the other input nodes to the highly conductive pathways.This behaviour was also found in hardware experiments (e.g.Ewers, 1982;Ford et al., 2000).Finally, inflow entering at the recharge side is routed to the discharge side almost completely through large-aperture conduits.Thus, main pathways continue to enlarge, whereas the remaining protoconduits do not evolve much further.
this late stage the fraction of tubes with turbulent flow decreases with ongoing conduit enlargement.This is because the average flow velocity v decreases faster than the hydraulically effective aperture a increases if the discharge stays constant ( ).According to Eq. ( 1) this leads to lower Reynolds numbers.
The general evolution and the resulting conduit patterns correspond well to those from earlier investigations that were based on similar or the same dissolution rate laws (e.g.Lauritzen et al., 1992;Howard and Groves, 1995) but also to those of Bloomfield et al. (2005), who used a different, generic type of dissolution rate law.This also applies to the effects of the intial aperture variability.For instance, lower initial heterogeneity results in lower tortuosity of the evolving pathways as compared to the scenario shown in Fig. 2, which agrees with findings by Howard and Groves (1995) and Bloomfield et al. (2005).Thus, initial heterogeneity provides an important control on the resulting conduit pattern.bimodal aperture distributions evolve at mature stages after breakthrough.The mean aperture sizes of strongly enlarged tubes in the high-heterogeneity models are slightly smaller than those of the low-heterogeneity models, which is probably an effect of the different breakthrough times (see below).Nevertheless, the effect of initial protoconduit heterogeneity on mean aperture sizes of mature conduits is minor.However, the number of largely widened tubes strongly depends on initial heterogeneity: their quantity decreases with increasing initial aperture variability.This confirms findings by Bloomfield et al. (2005), which were obtained with a generic dissolution rate law instead of the carbonate dissolution kinetics considered here.

Aperture distributions
The temporal evolution of the mean aperture distribution and that of the individual realisations for the highheterogeneity models is shown in Animation 2. Before breakthrough there is a slight widening of all apertures and the distribution stays unimodal for all realisations.Immediately after breakthrough some pathways connecting the recharge and discharge sides attract almost all the water (Animation 1), widen quickly and thus create a second peak in the aperture distribution, whereas the widening of nearly all tubes outside these preferred flowpaths ceases.
Breakthrough times range from approximately 5 ka to more than 30 ka in the high-heterogeneity settings but only from approximately 2.5 ka to 5 ka in cases of low initial protoconduit heterogeneity.Thus, increasing initial heterogeneity leads to higher variability of predicted breakthrough times.In addition, the mean breakthrough time in the high-heterogeneity settings is evidently higher than that in the low-heterogeneity settings.This appears to contradict results by Hanna and Rajaram (1998), who found lower breakthrough times with increasing initial heterogeneity.
These authors, however, attempted to compare settings with identical initial bulk permeability, whereas in the simulations reported here the mean of the logarithm of the apertures is kept constant.As illustrated by Fig. 3 there are much more small apertures in the initial aperture distribution of the high-heterogeneity settings than in the low-heterogeneity settings.For instance, the median of the highly heterogeneous aperture distribution is approximately 0.35 mm (average of all realisations), whereas it is 0.5 mm in the low-heterogeneity settings.The hydraulic conductivity of the protoconduits decreases strongly with decreasing aperture and thus there are more low-conductive protoconduits in the heterogeneous setting.Since the initial bulk permeability is controlled by the low-conductive protoconduits, it tends to decrease with increasing heterogeneity, which explains the observed differences in the breakthrough times.
It is further noteworthy that initial flow rates vary strongly among the 30 realisations of the heterogeneous setting (cf.Animation 2).This illustrates that protoconduit networks with (nearly) identical aperture distributions may have very different bulk permeabilities depending on the spatial arrangement of the individual apertures, which finally lead to different breakthrough times.Despite the high variablity of breakthrough times in the realisations of the heterogeneous model, the evolving large-aperture conduits approach similar sizes after breakthrough (cf.Animation 2).A similar finding is reported by Birk et al. (2000) for single conduits that are enlarged by dissolution following a diffusion-controlled firstorder rate law: while the breakthrough time of the conduit was found to be highly sensitive to the initial aperture, the long-term development was insensitive to this parameter.Thus, in the long term the initial heterogeneity of the protoconduits controls the number of large-aperture conduits that evolve but it hardly affects the size of these conduits.
In accordance with the finding by Birk et al. (2005), the evolution of the conduit pattern shown in Animation 1 suggests that the structure of a conduit network is determined at an early stage before breakthrough.At the mature stage after breakthrough the large-aperture pathways attract most of the water and continue to grow, while the other tubes remain nearly unchanged.These two evolutionary stages and the quick separation of continuously growing large-aperture conduits from non-evolving apertures after breakthrough can also be identified in the evolution of the aperture distribution: the breakthrough labels a splitting of the initially log-normally distributed apertures into two separated distribution modes, i.e. a distinct gap ranging over some orders of magnitude separates the fractions of smallsized and strongly enlarged tubes.The resulting bimodal aperture distribution is dynamically stable and evolves both in the low and high heterogeneous settings.Similar findings were reported by Bloomfield et al. (2005) and Rehrl et al. (2008Rehrl et al. ( , 2010) ) considering a generic dissolution rate law and gypsum dissolution kinetics, respectively.

Influence of recharge limitations
This section considers the influence of the maximum available recharge on the evolving aperture distribution.To this end, Q max is varied by steps of one order of magnitude.
For values of Q max ≤ 10 −9 m 3 s −1 , conduit enlargement is negligible within a simulation period of 30 ka.For Q max ≥ 10 0 m 3 s −1 , some tubes exceed aperture widths of 10 m, and thus parallel tubes potentially will be merged, which is not considered by the model employed here.Therefore, scenarios with maximum recharge rates ranging from 10 −8 m 3 s −1 to 10 −1 m 3 s −1 are analysed.
Figure 4 shows the evolving mean aperture distributions in settings with identical (low) degree of initial heterogeneity but four different maximum recharge rates.
The corresponding temporal development of mean aperture distributions at all eight maximum recharge rates is illustrated in Animation 3.
Dependent on the amount of available water, the evolving aperture distributions suggest a distinction between two different types of conduit development, which are schematically sketched in Fig. 5.If the total flow rate can increase to sufficiently high values, breakthrough occurs.The development of conduit patterns in these scenarios is similar to that described in the previous section.At the initial stage the protoconduits compete for flow, and those with large apertures and high hydraulic gradients (or more generally with low breakthrough times) succeed in attracting the majority of flow.After breakthrough they widen under fast first-order kinetics and separate from the bulk distribution, thus creating bimodal aperture distributions.In the following, this type of conduit development is termed competitive flow.If the maximum recharge is restricted to some values below a critical flow rate (≈Q bt in Fig. 5, which is ≈10 −6 m 3 s −1 in the scenarios presented here), no clear separation of large-aperture conduits from a fraction of non-evolving protoconduits takes place.The aforementioned competition for flow is suppressed and the positive feedback mechanism between increasing flow rate and dissolutional widening is switched off, which is why we term this conduit development under suppressed flow conditions.As suggested by Fig. 5, intermediate cases in between these two types of conduit development have to be considered, too.

Competitive flow
This section refers to scenarios in which conduits evolve under evidently competitive flow conditions leading to dynamically stable bimodal aperture distributions, i.e. for simulations with Q max ≥10 −6 m 3 s −1 .An aperture size of approximately one centimetre marks the threshold between shown.If Q in is restricted to ≈Q bt or higher, competitive flow systems establish and bimodal aperture distributions develop.In suppressed flow systems the total flow rate is limited to some lower value and lead to unimodal aperture distributions at mature stages.Intermediate cases must be considered.
tubes which stay small-sized and those which continue to grow.This threshold applies to both high-heterogeneity (cf.Animation 2 and Fig. 3) and low-heterogeneity (cf.Animation 3 and Figs. 3, 4, and 6) settings.
As opposed to the initial heterogeneity, the maximum available recharge is found to have an effect on the aperture size of the largest conduits evolving in these scenarios (Fig. 4): the first-order dissolution rate in the pathway that has broken through is proportional to e − 1 Q and the aperture growth rate is proportional to the dissolution rate (e.g.Dreybrodt et al., 2005).Hence, the aperture growth rate increases with increasing flow rate and consequently the higher the maximum flow rate the larger the maximum aperture size.The limitation of recharge also has a significant effect on the number of tubes evolving: the number of tubes that do not evolve increases with decreasing maximum recharge.Thus, the higher the available recharge the more large-aperture tubes evolve.
Snap-shots of typical mean aperture distributions at different stages of competitive conduit development are shown in Fig. 6.
Animation 4 shows the temporal development of the corresponding mean distributions as well as that of the individual realisations.stage of conduit development before breakthrough.After breakthrough, which is marked by the separation of a fraction of evolving large-aperture tubes, the distribution of smallaperture tubes remains nearly unchanged and still appears to be reasonably well described by a log-normal distribution (Fig. 6 and Animation 3).The total number of evolving large-aperture tubes stays constant during the mature stage of conduit development.This follows from the fact that the distributions of small-aperture tubes remain unchanged after breakthrough.The evolution of the aperture distributions in the lowheterogeneity realisations shown in Animation 4 further reveals that the random intial heterogeneity causes slightly different breakthrough times and thus differences in the aperture distributions evolving early after breakthrough.Nevertheless, the fractions of strongly widened tubes converge to similar distributions quite soon for all realisations, although the appearance of the individual conduit pattern is very different.Thus, the resulting long-term aperture distribution is determined by statistic properties of the initial aperture distribution (here, µ 0 and σ 0 of the initial log-normal distribution) under the given hydrogeological conditions (e.g.maximum recharge rate) rather than by the specific arrangement of the protoconduits within the individual networks.This is further confirmed by Fig. 7, which compares for all realisations of selected scenarios with high initial heterogeneity the percentage of large-aperture tubes after 30 ka with the breakthrough times and the time when the based on breakthrough based on inflow limitation Q max = 10 −4 m 3 s −1 Q max = 10 −5 m 3 s −1 Q max = 10 −6 m 3 s −1 Q max = 10 −7 m 3 s −1 Fig. 7. Comparison of the percentage of large-aperture tubes (a≥0.01 m) after 30 ka with breakthrough time (filled symbols) and time when the constant maximum inflow rate is reached (open symbols) in selected high-heterogeneity scenarios ( σ 0 µ 0 = 1).Note that for Q max = 10 −4 m 3 s −1 and Q max = 10 −5 m 3 s −1 both times practically coincide and that in the scenario with Q max = 10 −6 m 3 s −1 breakthrough occurs only in some realisations.Breakthrough does not occur in any realisation with Q max = 10 −7 m 3 s −1 .constant maximum flow rate is reached, respectively.While the percentage of strongly widened tubes is clearly dependent on Q max , no correlation with the breakthrough times or the times when the maximum flow rate is reached under competitive flow conditions is observed.The figure further reveals that with increasing Q max the percentage of the largeaperture tubes varies more strongly between the individual realisations.Hence, the more recharge is available the more important becomes the specific arrangement of the protoconduits.In contrast, if Q max is sufficiently low the resulting percentage of strongly widened tubes varies only within a narrow range and thus is mainly determined by the available amount of recharge.The time when the maximum flow rate is reached is nearly identical to the breakthrough time in all scenarios where breakthrough occurs.However, in most cases of Q max = 10 −6 m 3 s −1 breakthrough does not occur, i.e. calcium concentrations exceed the switch concentration before the water reaches the discharge side of the model domain.Nevertheless, bimodal aperture distributions are found to evolve in all realisations of the scenario (Animation 4), which suggests that competitive conduit development does not necessarily involve breakthrough.If the maximum available recharge is further reduced breakthrough does not occur in any realisation.These scenarios are presented in the following section.

Suppressed flow
Animation 3 and Figs. 4 and 7 suggest that the evolution of aperture distributions in the two scenarios with very strongly limited recharge rates, i.e.Q max = 10 −7 m 3 s −1 and Q max = 10 −8 m 3 s −1 , differs from that described in the previous section.As opposed to the scenarios with higher maximum available recharge, no fraction of large apertures clearly separates from the bulk of small-aperture tubes.The distance at which the switch to the slow higherorder kinetics is reached increases linearly with the flow rate.In the scenarios considered here, recharge is so strongly limited that the increase of flow rates is suppressed at a level that does not permit the calcium concentration at the discharge boundary to fall below the switch concentration.Hence, the positive feedback between conduit enlargement and increasing flow rates is switched off before breakthrough and thus breakthrough does not occur.As a consequence, under these suppressed flow conditions conduit development is found to be non-competitive, i.e. all tubes continue to enlarge during the entire simulation period and the aperture distribution remains unimodal rather than bimodal.The difference between conduit development under competitive flow and suppressed flow conditions is most evident in the different evolution of the aperture distributions in the time period after breakthrough (see Animation 3).The distribution of the small apertures remains unchanged after this time under competitive flow conditions (see previous section) but continues to change in the two scenarios considered here.
In order to assess whether or not bimodal aperture distributions evolve in the long-term even at strongly limited recharge rates, simulations extending to 3 Ma were performed for Q max values of 10 −7 m 3 s −1 and 10 −8 m 3 s −1 and both degrees of initial heterogeneity.Animation 5 shows the long-term temporal development of aperture distributions of all realisations with low initial heterogeneity for both maximum recharge rates.The mean aperture distribution of the realisations with low initial heterogeneity and Q max = 10 −8 m 3 s −1 at selected simulation times are shown in Fig. 8.One example of corresponding conduit pattern at a late stage is given in Fig. 9.
In the long-term simulations with Q max = 10 −8 m 3 s −1 , the calcium concentration is very high in all tubes and exceeds 0.99c eq already at the outlet of the tubes connected to the recharge boundary.Thus, dissolution rates are very low throughout the whole network and the entire range of apertures continues to be slowly enlarged.Hence, the aperture distributions stay unimodal even after time periods of 3 Ma both in the case of low (Fig. 8  runs suggest that the distribution does not become bimodal even after several hundred million years.Instead, the distribution becomes uniform.Nevertheless, in Fig. 8 a small fraction of large tubes is slightly separated from the bulk distribution.However, this separation occurs at a very early stage (t < 10 ka) and subsequently the large apertures do not separate further from the bulk distribution.Indeed, the largest tubes are directly connected to the recharge boundary and thus receive inflow with zero calcium concentration.Because of the very low flow rates, the calcium concentration approaches values close to the equilibrium concentration c eq still within these tubes (cf.Fig. 9), and thus the dissolution switches from first-order to higher-order kinetics before reaching the outlet.The operation of first-order kinetics close to the inlet yet is sufficient to cause an initially faster growth and thus the slight separation in the aperture distribution.This effect appears in the other scenarios, too (cf.Figs. 3  and 4, and Animations 2 and 3).
As opposed to the above scenario with very strongly limited inflow, the long-term simulations with Q max = 10 −7 m 3 s −1 show a behaviour in between the two endmembers obtained under competitive and suppressed flow conditions.
Thus, the scenario corresponds to the intermediate case indicated in Fig. 5.An example of conduit pattern after 3 Ma is illustrated in Fig. 10.Because of the low flow rates, at the initial stage all tubes are widened slowly.Therefore, at the early stage the entire range of apertures continues to enlarge evenly and the distribution stays unimodal (cf.Animation 5).After some 3.00Ma time a few channels have invaded deeply into the aquifer but without reaching the output boundary.Tubes along these channels exhibit slightly lower concentrations and thus higher dissolution rates.As a consequence, preferred pathways growing slowly towards the discharge boundary develop (Fig. 10).This becomes obvious in the aperture distribution after several ten thousand years (the time varies between the different realisations), which in the range of millimetres and below becomes stable while larger tubes still grow.From then on the aperture distribution exhibits a minimum at aperture sizes of approximately 1 cm.Only tubes that exceed this size, i.e. those forming the preferred pathways, grow further and therefore the number of tubes with apertures in the order of 1 cm decreases.Thus, the aperture distribution changes towards a bimodal distribution just as in the simulations with higher flow rates.However, the resulting distribution is different from those obtained with higher flow rates in that even the initially smallest apertures have been affected by considerable dissolutional widening.It is further noteworthy that the maximum of the aperture distribution obtained with Q max = 10 −7 m 3 s −1 is found to be at a lower aperture width than that obtained with Q max = 10 −8 m 3 s −1 .Thus, the higher flow rate causes the development of a few prominent pathways but the apertures of the bulk distribution are larger with the lower flow rate (cf.Animation 5 and Figs. 9 and 10).

Discussion
Figure 11 provides an overview of the temporal evolution of the percentage of large-aperture tubes in simulations with different degrees of initial heterogeneity and different maximum recharge rates.Corresponding statistical values are given in Table 2.
The simulation results confirm that the general findings by Bloomfield et al. (2005), who assumed a generic growth law in their simulations, are also valid if conduit enlargement is governed by the empirically established rate laws for carbonate dissolution.In particular, it is found that bimodal aperture distributions with a constant percentage of largeaperture tubes evolve after breakthrough if the recharge is not strongly limited.Only these large-aperture tubes continue to grow, while the apertures of the other tubes Hydrol.Earth Syst.Sci., 15, 3715-3729, 2011 www.hydrol-earth-syst-sci.net/15/3715/2011/ 3.00Ma remain nearly unchanged.In agreement with Bloomfield et al. (2005) we find that the number of evolving large-aperture tubes decreases with increasing heterogeneity of the initial aperture distribution (Fig. 11).Following the terminology introduced by Worthington (2009) settings with high initial heterogeneity thus favour macrokarst, those with low initial heterogeneity microkarst.As outlined in the introduction, numerical simulations of the evolution of conduit pattern have revealed additional factors that favour the development of macrokarst, such as steep hydraulic gradients or low chemical saturation of the recharge.It is thus likely that the aperture distribution is influenced by these factors in a way similar to the effect of initial heterogeneity described above.
In this work one additional factor is examined that has received little attention before, namely the limitation of the recharge to the soluble unit.Our results demonstrate that the number of large-aperture tubes that continue to grow after breaktrough decreases with decreasing maximum recharge (Fig. 11), which corresponds to Kaufmann and Braun (1999) who found that in model scenarios with fixed recharge more branches develop the higher the recharge rates.In this regard, the effect of limiting the maximum recharge is similar to that of increasing heterogeneity of the initial aperture distribution.However, the effect on the variability of the number of large-aperture tubes evolving in the realisations of the individual scenarios is different (Fig. 11): the variability of the resulting percentage of large-aperture tubes is found to decrease with decreasing available recharge but nearly unaffected by the degree of initial heterogeneity.In addition, the differences between the scenarios with low and high degrees of initial heterogeneity reduce with decreasing available recharge.Hence, the lower the amount of available recharge the more it controls the number of large-aperture tubes that evolve.
The simulation results summarized in Fig. 11 further reveal that the type of conduit development changes more fundamentally if the recharge is very strongly limited.The feedback between conduit enlargement and increasing flow rate, which is operating under constant-head boundary conditions, is then found to be suppressed due to the limited availability of recharge.As a consequence, there is a slow, ongoing enlargement of all pathways such that the percentage of large-aperture tubes continues to increase during the entire simulation time.As opposed to the  bimodal aperture distribution evolving under competitive flow conditions, the aperture distribution is found to stay unimodal if the increase of flow is suppressed at very low recharge rates.
Although the availability of water is inherently limited in any type of hydrogeologic setting, the question arises if and where the recharge limitation is strong enough to create unimodal aperture distribution under suppressed flow conditions.In unconfined settings, the maximum recharge is determined by precipitation and evapotranspiration.Yet the network of interconnected conduits that evolves in the subsurface increases with time and thus the size of the catchment area increases.Although the maximum size of the catchment is limited by hydrogeologic boundaries such as impermeable rocks the amount of available water in these systems generally is expected to be sufficient for competitive conduit development and thus the evolution of bimodal aperture distributions.In contrast, confined flow in artesian basins is expected to be more restricted than flow in unconfined settings.Simulations by Rehrl et al. (2008Rehrl et al. ( , 2010) ) addressing the evolution of conduits in artesian gypsum settings of the Western Ukraine indeed suggest that the flow rates may be sufficiently low to cause the development of unimodal aperture distributions.Because of the high solubility of gypsum, a few hundred thousand years were found to be sufficient to create maximum aperture sizes in the range of metres, whereas the maximum apertures occurring in the carbonate setting considered here are in the centimetre range even after 3 Ma (Fig. 8 and Animation 5).Nevertheless, conduit development under suppressed flow conditions may constitute an important mechanism in the hypogenic speleogenesis in deep, confined carbonate aquifers that causes significant modifications of the hydraulic aquifer properties.The occurrence of microkarstic aperture distributions such as those shown in Fig. 8, for instance, may be an important aspect in the assessment of thermal resources (cf.Goldscheider et al., 2010) or the transport and entrapment of hydrocarbons in deep carbonate aquifers (Klimchouk, 2007).
Table 2. Statistical values of the percentage of large-aperture tubes with a≥0.01 m after fixed genesis times for different maximum recharge rates Q max (m 3 s −1 ) and both degrees of initial heterogeneity σ 0 µ 0 .N gives the number of realisations from which the mean percentage µ and the corresponding standard deviation σ are calculated.Realisations in which nary a tube reached the given size range are excluded.
after 30 ka: Suppressed conduit development in deep, confined aquifers may also constitute an important precondition for conduit development under unconfined conditions after uplift and erosion of the confining layers.The general characteristics of suppressed conduit development identified in our model simulations, namely slow and uniform enlargement of the entire set of pathways during large time periods, fit well to those associated with the inception horizon hypothesis advanced by Lowe (1992Lowe ( , 2000)).The inception horizon hypothesis assumes that the earliest stage of conduit development is related to horizons within a carbonate sequence that are particularly prone to dissolution.The occurrence of inception horizons that favour the development of solution conduits was recently supported by field evidence (Filipponi, 2009).Within the framework of the inception horizon hypothesis, conduit evolution under suppressed flow conditions (cf.Fig. 8) may represent the inception phase, which is later (e.g. after uplift and erosion) followed by a gestation phase (Lowe, 2000) with competitive conduit development.Thus, suppressed flow and competitive flow conduit development might represent successive stages of conduit evolution rather than independent types of evolution in different hydrogeological settings.
The simulations in this work were based on the empirically established rate laws for carbonate dissolution.Yet the results referring to competitive conduit development were found to be in good agreement with published results from simulations that were based on different growth laws.Thus, we are confident that even in the case of suppressed flow conditions the general conclusions are not strongly affected by the type of dissolution rate law.Nevertheless, the influence of the various dissolution mechanisms and sources of the chemical aggressiveness of the water deserve further consideration, in particular, with respect to deep carbonate aquifers.For instance, dissolution can be enhanced by mixing corrosion (Dreybrodt and Gabrovšek, 2000) or the occurrence of geogenic acids (e.g.Goldscheider et al., 2010) and thus conduit development under suppressed flow conditions might be more efficient than in the simulations reported here.

Conclusions
The numerical simulations presented here reveal that the dissolutional enlargement of discrete pathways in carbonate rock generally creates bimodal aperture distributions where only the fraction of large-aperture pathways continues to grow after the breakthrough of a prominent passage.The number of large-aperture conduits that continue to grow is found to decrease with increasing heterogeneity of the initial aperture distribution and with decreasing availability of recharge.The type of conduit development is found to change fundamentally if the maximum available recharge is severely reduced.Under these suppressed flow conditions unimodal aperture distribution evolve and the entire range of apertures continues to grow slowly, resulting in rather uniform aperture distributions in the long-term.Although this type of conduit development is not efficient in creating large voids, it does succeed in creating a multitude of pathways with apertures in the order of millimetres up to centimetres within millions of years.Since this appears to be sufficient to change the hydraulic properties of the rock significantly, it needs to be considered in the hydrogeological assessment of deep carbonate aquifers.
Our simulation results further show that the breakthrough time and the mature conduit pattern may vary strongly between realisations with identical statistical properties, whereas the aperture distributions evolving after breakthrough are generally found to be very similar.In particular, the dependency of the resulting aperture distribution on the initial variability of aperture sizes and their spatial arrangement reduces with decreasing availability of recharge.Thus, the aperture distribution of mature karst systems is found to be determined by hydrogeological factors such as initial heterogeneity and availability of recharge but rather insensitive to the random spatial arrangement of individual apertures within the network of protoconduits.Hence, the analysis of aperture distributions observed at a given site might be well suited to provide additional insight into the environmental conditions that controlled the speleogenetic evolution.

Fig. 1 .
Fig. 1.Conceptual model of a conduit network, potentionally developing along a bedding plane embedded in fractured soluble rock.The modelling domain represents a two dimensional section within a limestone aquifer.Recharge into and discharge from the system is restricted to two opposite sides, whereas the remaining sides are no flow boundaries.

BFig. 2 .
Fig. 2.Example of a conduit pattern at a late stage (30 ka) of karst evolution after breakthrough showing one realisation with high initial protoconduit aperture heterogeneity and an intermediate total maximum recharge rate of Q max = 10 −3 m 3 s −1 .Filled circles represent nodes.Rings mark nodes located at the recharge side which reached maximum inflow rates of q max = 5×10 −5 m 3 s −1 .Single tubes are shown as arrows indicating flow direction.The width of the arrows represents the aperture size and tildes denote regions with turbulent flow.Red colour indicates regions where first-order dissolution kinetics is active.

Figure 3
Figure3compares mean aperture distributions after 30 ka of karstification from 30 realisations with identical (intermediate) maximum flow rate but low and high initial aperture heterogeneity, respectively.In either case

Fig. 4 .
Fig. 4. Mean aperture distributions at late stages of karstification in models with low initial aperture heterogeneity and different maximum recharge.

Fig. 5 .
Fig. 5. Schematic development of the total inflow rate Q in in time.The time (t bt ) and flow rate (Q bt ) at breakthrough (bt) areshown.If Q in is restricted to ≈Q bt or higher, competitive flow systems establish and bimodal aperture distributions develop.In suppressed flow systems the total flow rate is limited to some lower value and lead to unimodal aperture distributions at mature stages.Intermediate cases must be considered.

Fig. 6 .
Fig. 6.Temporal evolution of the mean aperture distribution of the low-heterogeneity model with total maximum recharge Q max = 10 −4 m 3 s −1 .

Fig. 8 .
Fig. 8.Long-term evolution of the mean aperture distribution in the low-heterogeneity settings with strongly limited maximum recharge rates of Q max = 10 −8 m 3 s −1

B
. Hubinger and S. Birk: Influence of hydrogeological conditions on aperture distributions

Fig. 9 .
Fig. 9. Example of conduit pattern showing the same realisation illustrated in Fig. 2 but for the long-term simulation of suppressed flow scenarios with Q max = 10 −8 m 3 s −1 and low initial aperture heterogeneity.

Fig. 10 .
Fig. 10.Conduit pattern of the realisation shown in Fig. 9 but for the intermediate case with Q max = 10 −7 m 3 s −1 .

B
. Hubinger and S. Birk: Influence of hydrogeological conditions on aperture distributions

Fig. 11 .
Fig. 11.Overview of the percentage of large-aperture tubes with a ≥ 0.01 m evolving in simulations with different initial heterogeneity and different maximum recharge rates.The left panel (a) shows the temporal evolution in a network with low initial heterogeneity for one realisation with Q max ≥10 −6 m 3 s −1 (until 30 ka) and for all realisations with Q max ≤10 −7 m 3 s −1 (suppressed flow scenarios; until 3 Ma), respectively.Dashed lines represent the supposed further development of the individual realisations until 10 Ma of karstification.The panels (b) and (c) illustrate results for both degrees of intial heterogeneity after 30 ka excluding realisations in which nary a tube reached an aperture width of a≥0.01 m.The corresponding long-term results for suppressed flow scenarios after 3 Ma are marked with arrows.

Table 1 .
Parameter values used for the simulations.