D. Romanov ^{(1)}, F. Gabrovsek ^{(2)} and W. Dreybrodt^{ (3)}

^{(1)} Karst Processes Research Group, Institute of Experimental Physics, University of Bremen,

28359 Bremen, Germany

^{(2)} Karst Research Institute ZRC SAZU, Postojna, Slovenia

^{(3)} Karst Processes Research Group, Institute of Experimental Physics, University of Bremen, 28359 Bremen, Germany

## Abstract

First models of karst evolution considered a single isolated fracture with no loss of flow along its entire length. Under conditions of constant head dissolution of limestone creates a positive feedback-loop of increase of aperture widths and flow until at breakthrough the flow and aperture width are enhanced dramatically. If a second dimension is added to this model domain, in the simplest case by an exit-tube connected to the isolated channel, water loss from the isolated channel occurs. We have investigated the influence of the water loss on the breakthrough time of the single channel. In all cases, when water loss is present, more aggressive solution enters at the input. The aggressive solutional activity penetrates deeper along the conduit. Therefore dissolutional widening at the exit is enhanced and breakthrough times are reduced. This is discussed in detail by investigating the profiles of hydraulic head, flow rates, aperture widths, and calcium concentrations along the conduit as they evolve in time and comparing them to those of the isolated 1-D conduit.

In a further step the 1-D conduit is embedded into a net of fractures with smaller aperture widths. The conduit is located in the center of the rectangular domain and connected to the 2-D net at equally spaced nodes. By this way exchange flow from the conduit into the net can arise. But also flow from the net to the conduit is possible. We have studied the evolution of this aquifer considering dissolution also in the network of the narrow fissures. Flow from the main central fracture into the net again reduces breakthrough times. After breakthrough, however, a complex exit fan evolves in the net, which later on is overprinted by a net of entrance fans propagating down flow. These fans are related to flow from the net into the central fracture. The evolution of these fans resulting finally in a maze-like structure is significant for high hydraulic gradients (i ³ 0.1) as they exist at artificial dam sites. For such situations realistic modeling has to include dissolutional widening in the net. For low hydraulic gradients, i<0.03, the evolution in the net is slow compared to that of the central conduit and therefore the aquifer is dominated by the evolution of the central fracture.

Keywords: Limestone, Karst aquifers, Speleogenesis, Modeling, Conduit evolution.

## 1. Introduction

The evolution of karst terrains is governed by many factors, such as climate, geological settings, the location and geometry of the catchment areas, the type of the soluble rock, and the distribution of primary fractures in this rock. Karst has been a subject of extensive research since centuries. The articles of Shaw (2000), Lowe (2000), and White (2000) are interesting reviews about the development of the speleogenetic studies from ancient times to present days. People are interested in karst evolution not only because of the beauty of the karst landforms, but also because of their practical importance. Karst aquifers (rock bodies sufficiently permeable to transmit groundwater (Bear and Veruijt, 1987) are the main source of drinking water for about 25% of the world population (Ford and Williams, 1989). Sinkholes, sinking streams, closed depressions and caves characterize the topography of karst terrains. All these different landforms have a common element – the well developed subsurface drainage system.

Initially, when the hydraulic conductivity of the soluble rock is small, only a small amount of water can enter. The initial aperture widths of the fractures are in the range of several 100 m m . Therefore the flow through them is laminar. CO2 containing water is an aggressive solution capable to dissolve limestone. When it flows through the fissures in the soluble rock their aperture widths increase due to dissolutional widening. Some of the fissures widen faster than others. The flow through these increases and consequently the rate of their widening becomes higher. This positive feedback loop is the reason for the development of secondary porosity and consequently for the development of a complex, extremely heterogeneous aquifer. Flow through some of the widened fractures finally becomes turbulent. Until this time the hydraulic conductivity of the karst aquifer has increased by orders of magnitude.

Nowadays, we have the basic knowledge about the properties of karst aquifers, about the hydrodynamics of flow through them, and about the dissolution kinetics of the soluble rocks. Together with the computational power, this enables us to build numerical models of karst aquifers and study them. Dreybrodt (1988, 1990, 1996) and Palmer (1988, 1991) have presented the first numerical models. They are constructed on the basic principles of groundwater chemistry and hydrology, and study the evolution of an isolated one-dimensional conduit under various boundary conditions.

Using the information about the evolution of a single fracture, we are able to build and understand more complex two-dimensional models. Lauritzen et all. (1992), Groves and Howard (1994), and Howard and Groves (1995) presented models for the evolution of two-dimensional networks. Siemers and Dreybrodt (1998), Siemers (1998), and Dreybrodt and Siemers (2000) have investigated the evolution of two dimensional percolation networks under various lithological and hydraulic conditions. They extended their studies for cases of practical interest, namely the karstification in the vicinity of large hydraulic structures. The karstification below dam sites is investigated by Dreybrodt et al. (2002) and Romanov et al (2003). Clemens et al (1997a; 1997b; 1996), and Bauer (2002) reported on a double porosity model. They couple the large conduit flow with the flow in the surrounding continuum of narrow fissures and calculate the evolution of the conduits. Kaufmann and Brown (1999, 2000) and Kaufmann (2003) used a similar approach. In their model, however, prominent conduits are incorporated into the continuum.

Gabrovsek (2000), Gabrovsek and Dreybrodt (2000a, 2000b) have studied the evolution of a single fracture and two dimensional percolation networks under various chemical and hydrological boundary conditions. Together with the numerical results, several analytical estimations for the breakthrough time are given. Gabrovsek and Dreybrodt (2001) have put forward also a model for the evolution of an unconfined aquifer.

Any of these different modeling approaches has its advantages and disadvantages. An important result is that all of them give similar results for basic scenarios, specially designed for comparison.

The goal of this paper is to discuss the transition from one-dimensional (single fracture) models to two-dimensional (fracture networks) ones, and the new mechanisms arising. Therefore we start with a short overview of the evolution of a single fracture under a constant head boundary condition. We compare these results with the evolution of the same fracture, but embedded into a two-dimensional network.

## 2. Theoretical background

The basic element of our models is the single conduit. We use this block to create complex systems and model the processes, which govern the evolution of natural karst aquifers.

A single conduit is presented on Fig.1. It has a rectangular shape, but it could be a cylindrical tube or have any kind of characteristic geometry. The aperture width , the width , and the length *L*, characterize the rectangular fractures.

** Fig. 1.** Single rectangular conduit, characterized by its length *L*, width , and aperture width .

It is possible also to model the evolution of cylindrical conduits (tubes), but for the rest of this work we will deal only with rectangular fractures.

To model the evolution of the fracture in time, we need to know:

- The hydrological laws – governing the flow through the fracture;
- The chemical laws – governing the change of the calcium concentration along the fracture caused by dissolution.

Details about these elements are given in the literature: (Dreybrodt, 1988; Dreybrodt, 1996; Siemers and Dreybrodt, 1998; Gabrovsek, 2000).

** Fig. 2.** Mass conservation in the part of the fracture between *x* and

To calculate the widening rate *F(x)* we need to know the concentration *c(x)* along the fracture. Fig. 2 represents part of a fracture between *x* and *x+** D x* , where *x* is the distance from the entrance. We use conservation of mass and obtain the following equation:

(1) ,

where *A(x)* is the cross-sectional area at *x*, *P(x)* the perimeter *v* is the velocity of the fluid and *Q* is the flow rate. The solution of this equation for plane parallel walls is given by (Dreybrodt, 1996; Dreybrodt and Gabrovsek, 2000; Gabrovsek, 2000) (2):

where *c in* is the concentration at the input of the fracture and:

(3)

* P* is the perimeter, *k 1* and *k n* are kinetic constants, *c eq* is the equilibrium concentration with respect to calcite, *n* is the kinetic order, and *x s* is the position of the switch between the linear and the non-linear rate law (see Eisenlohr et all. 1999). Note the importance of the nonlinear part for the karstification. If the rate law is fully linear then the rates decrease exponentially along the fracture. By this way the exit of the fracture remains practically unaffected. On the other hand, because of the non-linear rate law, dissolution is active along the whole length of the fracture. By this way its exit part is widened, and flow increases. Consequently the dissolution rates increase, the exit is widened faster, and the flow rises faster in time. This positive feedback loop leads to the breakthrough event. The moment of this event is called - breakthrough time. The flow increases by several orders of magnitude at this moment, and the concentration at the exit becomes practically equal to *c in*. After breakthrough the fracture continues to widen evenly along its entire length.

It is not possible to give an analytical expression for the evolution of the fracture width with time. But if we assume that the rates at the exit are active along the entire length of the fracture, then its walls remain parallel during the entire evolution, a reasonable approximation by an analytical solution can be obtained. The widening along the conduit, is given by the equation:

(4),

where (5):

*a 0* is the initial aperture width, g converts the dissolution rates from mol × cm -2 × s -1 to retreat of bedrock in cm × year -1. g is 1.17 × 10 9 for limestone.

For details about solving the equation see (Gabrovsek, 2000). The result is:

(6).

The breakthrough time *T B* is:

(7).

If we insert Eq. 5, we obtain an upper limit of the breakthrough time *T B*, and also its dependence on the basic parameters determining karstification (Dreybrodt, 1996; Dreybrodt and Gabrovsek, 2000; Gabrovsek, 2000) (8):

where * r* is the density of water, * h* its dynamic viscosity. *g* is the Earth’s acceleration. To obtain exact results one must take into account that the shape of the fractures shows a funnel like profile. Details on this numerical procedure are given by Gabrovsek (2000), Dreybrodt (1996).

**Fig. 3.** Two-dimensional rectangular network, and junction of four fractures

We use the single conduit to build a 2D network. Fig. 3 depicts an example of a rectangular two-dimensional fracture network. We are able to apply different hydrological and chemical parameters to every single fracture. Varying the spacing (*s*) between the fractures and their initial aperture widths we are able to model different hydraulic conductivities in the domain. The relation between the hydraulic conductivity and the parameters of the fractures is (Lee and Farmer, 1993):

(9),

where *s* is the spacing between the fissures.

A detailed description of the algorithms used to create 2D networks is given by Dreybrodt and Gabrovsek (2002).

## 3. Evolution of a single fracture

A single or isolated fracture is one, which cannot gain or loose flow along its length. The only way for it to exchange fluid with other fractures or a porous media is through its entrance and exit.

** Fig. 4. ** Pathways percolating through a block of soluble rock:

a) The block is impervious. Pathway 1 cannot exchange flow with Pathways 2 and 3;

b) The block carries a system of fine fractures (aperture widths *a 0*). The percolation pathway 1 (aperture width *A 0* is able to exchange flow with Pathways 2 and 3.

Fig. 4a depicts percolating pathways embedded into a continuous and impervious block of soluble rock. It is obvious, that every straight part of any of the pathways can be defined as a single conduit with an entrance, which is the exit of the preceding fracture. But in the context of the definition above, the whole pathway 1 can be defined as a single conduit. Its length is equal to the sum of the lengths of all straight fractures, which form it. For the pathways 2 and 3 this is not true. They are connected at the points A, B and C, and therefore they can exchange flow. This has a significant influence on their evolution.

Imagine now, that the block is not impervious, but has a structure of fine fractures - Fig. 4b. Of course any of them has its own width and length, but for simplicity we assume that they are equal. Now the percolating pathways are part of a network. What makes them different from the surrounding fine fractures is their initial aperture width. Pathway 1 cannot be accepted as an isolated conduit anymore. The question is: How can the fine network influence the evolution of the percolating pathways?

** Fig. 5.** Basic settings for the single fracture

a) Single fracture consisting of two isolated fractures connected at Node 1. (Fracture 0 – starting at Node 0 and ending at Node 1, and Fracture 1 – starting at Node 1 and ending at Node 2). The basic parameters are: *H=100m, h 2=0m, L=742.5m, a 0=a 1=0.03 cm*. The basic relations are shown in the figure;

b) Electrical circuit equivalent to the fracture system of Fig. 5a. The hydraulic resistances of the fractures are replaced with resistors. The potential difference *U* between the nodes represents the pressure head differences, and the current – *I*, the flow rate through the fractures. The basic relations are shown in the figure.

We start with the evolution of the single fracture shown on Fig. 5. The input Node is Node 0, and the output one is Node 2. We need one more node for the exchange of flow somewhere in fracture (later in this discussion). Therefore Node 1 is introduced. It splits the initial conduit into two parts with the same initial aperture widths *a 0=a 1=0.03 cm*. If no flow is exchanged through Node 1, the system remains unaffected. Therefore we threat these two fractures as one with a length equal to the sum of both lengths. The parameters and the boundary conditions for this model are listed in Table 1.

** Table 1 **

We introduce a parameter . Varying its value, we are able to move the exchange node along the fracture. The length of Fracture 0 is then . The length of Fracture 1 is . The total length *L=L 0+L 1* is kept constant and only the position of the Node 1 is varied.

The electrical equivalent of the hydrological set up from Fig. 5a is depicted in Fig. 5b. The resistors *R 0* and *R 1*have the same numerical value as the hydraulic resistances of Fractures 0 and 1. The hydraulic resistance of a fracture is given by (Beek and Muttzall, 1975):

(10),

* M* is a geometrical factor defined as

- ellipsoidal shape and

- rectangular shape.

** Fig. 6.** Evolution of the aperture widths, hydraulic head, concentration and flow rate for the scenario of a single fracture with length *L*.

a) Evolution of the aperture width. The behavior is typical for the single fracture – fast widening close to the entrance and much slower opening of the exit zone;

b) Evolution of the pressure head distribution. The highest gradients are at the tip of the opened zone;

c) Concentration along the fracture in units of the equilibrium concentration *c eq* . The concentration is low only in a relatively small area close to the entrance, and after the end of the widened part is close to saturation;

d) Evolution of the flow rate. This is a typical breakthrough curve of a single conduit. The flow rate is increasing slowly during most of the evolution time and then jumping with several orders of magnitude at breakthrough time.

The color code is the same for a), b) and c) and is depicted in Fig. 6c. BT=781 years

Fig. 6 depicts the numerical results for the evolution of the single fracture from Fig. 5. It shows the profiles of the aperture widths (Fig. 6a), the pressure head distribution (Fig. 6b), the concentration profiles (Fig. 6c), and the flow rates (Fig. 6d). The colors represent time, depicted in years. There is no flow out through Node 1. Therefore the flow rates through Fracture 0 and Fracture 1 are equal (Fig. 6d). They are limited by the constriction at the exit of Fracture 1 (see Fig. 6a). The pressure head at *t=0* is decreasing linearly along the conduit and then is changing with the change of the aperture width. The main head loss per length is at the tip of the widened zone (see Fig. 6b). The concentration profiles are also typical for the single fracture case (see Fig. 6c). Because of the slow growth of the exit part of Fracture 1, the flow rate remains low for most of the evolution time. The solution entering the fracture is attaining high Ca concentrations. Shortly before breakthrough, when the rate of widening of the exit is higher, the fresh solution starts to penetrate deeper into the channel. After breakthrough (at 781 years), the concentration becomes low along the entire length of the fracture (Fig. 6c) and therefore widening is even along the entire length.

Our next step is to make the one-dimensional model from Fig. 5 – two-dimensional. We use the system depicted in Fig. 7. There is a third fracture starting at Node 1. By this way flow is exchanged through it. The system Fracture 0 + Fracture 1 cannot be regarded as a single fracture anymore. Varying the position of Node 1 and the parameters of Fracture 2 we are able to study the effect on the evolution of the system from Fig. 5. In this paper we concentrate mainly on the changes in the hydrology. Therefore we use a configuration like the one shown on Fig. 7 – one input and two output fractures. Furthermore we simplify the problem using boundary conditions, which allow only flow out from the system (Fracture 0+Fracture 1) through Node 1. This is the most simple setup to study.

** Fig. 7.** Basic settings for the system of three fractures used to study the influence of exchange flow on its evolution.

a) Three single fractures (Fracture 0 – from Node 0 to Node 1, Fracture 1 – from Node 1 to Node 2, and Fracture 2 – from Node 1 to Node 3) are connected at Node 1. The basic parameters are: *H=100m, h 2=h 3=0m, a 0=a 1=0.03cm, L=742.5m, b 0=1m*. The different colors of the arrows represent the difference in the Ca concentrations of the inflowing and the outflowing solution;

b) Electrical circuit equivalent to the fracture system of Fig. 7a. The hydraulic resistances of the fractures are replaced with resistors. The potential difference *U* between the nodes represents the pressure head differences, and the current – *I* the flow rate through the fractures. The basic relations are shown on the figure.

The parameters of Fracture 0 and Fracture 1 are the same as those from Fig. 5a. The settings of the third fracture – Fracture 2, are varied. Therefore we need to introduce two more parameters. These are and . The first one is used to define by *m** × L* the length of Fracture 2 - . The second one is used to relate its initial aperture width *a 2* to *a 0* - . The boundary conditions are as follows:

- Constant head
*H=100 m*at the entrance node of fracture 0 (node 0); - Constant head
*h 2*=*h 3=0 m*at the exit nodes of fractures 1 and 2 (nodes 2, 3); - The initial aperture widths of fractures 0 and 1 are the same,
*a 0*=*a 1=0.03 cm*. The initial aperture width*a 2*of Fracture 2 is defined as , where . Dissolution is neglected along Fracture 2 and therefore the value of*a 2*remains unchanged until the end of the simulations –*a 2(t)=const*; - All fractures have the same width
*b 0=b 1=b 2*; - The lengths of the fractures are:

Fracture 0 – , where ,*0<k<1*; Fracture 1 - ;

Fracture 2 – , where . - The Ca concentration of the inflowing water is
*c in=0 mol/cm 3*. - The equilibrium concentration with respect to calcite is
*c eq=2e-6 mol/cm 3*.

The electrical circuit, equivalent to the setup of Fig. 7a is depicted in Fig. 7b. The differences between the initial situations of both scenarios – the isolated one from Fig. 5a and the exchange one (Fig. 7a) – are easy to see. A group of two parallel resistors (*R 1* and *R 2*) is connected with *R 0* at Node 1. The effective resistance connected to *R 0* is not equal to *R 1* anymore but is given by the equation:

, (11)

If the resistance *R 2* is infinitely large, then the circuit is equivalent to the one from Fig. 5b. ( ). This is also true for the hydraulic resistance of Fracture 2. If it has an infinite large value ( , see Eq. 10.), then the system represented by Fig. 7a is the same as the one of Fig. 5a, and behaves like a single fracture. But what happens in the case when *R 2* is finite? According to Eq. 11 will be less than the smaller of both resistances *R 1* and *R 2*. Consequently, the value of the effective resistance connected to *R 0* will be smaller. Actually this is the reason, which causes the differences in the evolution of both systems (Fig. 5a, and Fig. 7a.),: The value of the hydraulic head at Node 1 is different for both scenarios. In the two-dimensional scenario (Fig. 7) the hydraulic head at Node 1 is lower than in the isolated case. Therefore the hydraulic gradient along Fracture 0 is steeper and consequently flow is higher. According to Eq. 5 the dissolution rate at the end of the conduit (Node 1) will be higher, and the evolution will be faster (Eq. 8). At the same time because of the smaller gradient along Fracture 1 – its evolution will be delayed. But we are interested in the evolution time of the whole system (Fracture 0+Fracture 1). One can easily see that the situation at Node 1 will change with the evolution of Fracture 0 and Fracture 1. It is very difficult to obtain an analytical solution for the breakthrough time of the system Fracture 0 + Fracture 1. Therefore we will solve the problem numerically.

Varying the values of *f*, *m* and *k* we study the effect of the position of the exchange node, and the amount of the exchanged flow, on the evolution of the single fracture from Fig. 5.

3.1. Influence of exchange flow

We start with the variation of *f*. The values for *m* and *k* are fixed (*m*=*k*=1/2). The value of *f* is varied in the range from *f* almost equal to zero (this scenario is similar to the one presented by Fig. 5), to .

** Fig. 8.** Evolution of the profiles for a scenario with *f=0.33, k=1/2, m=1/2*.

a) Evolution of the aperture width. The behaviour is similar to the one of the isolated case (Fig. 5) – fast widening close to the entrance and much slower opening of the exit zone;

b) Evolution of the pressure head distribution. The highest gradients are at the tip of the opened zone;

c) Concentration along the fracture in units of the equilibrium concentration *c eq* . The concentration is low only in a relatively small area close to the entrance, and after the end of the widened part is close to saturation;

d) Evolution of the flow rate. Again a typical breakthrough curve of a single conduit is obtained. The flow rate increases slowly during most of the evolution time and then jumping with several orders of magnitude at breakthrough time. The flow through Fracture 0 and Fracture 1 is almost equal during the entire evolution.

The color code is the same for a), b) and c) and is depicted on Fig. 8c. BT=761 years

For *f*=0.33 Fig. 8 depicts the profiles for the aperture width (Fig. 8a), pressure head (Fig. 8b), concentration with respect to *c eq* (Fig 8c), and the evolution of the flow rate (Fig. 8d). The breakthrough time is only 20 years lower than the one for isolated case from Fig. 5. Profiles are very similar to those of the isolated fracture in Fig. 5. There is only one small detail, which is actually the reason for the reduction of the breakthrough time (see Fig. 8d). The flow rate through Fracture 1 is a little bit lower than the one through Fracture 0. The reason for this is the presence of Fracture 2. Its resistance is so high that it cannot significantly influence the evolution of the system. But even in this case we observe a reduction in the breakthrough time. We expect that this effect will become stronger with increasing values of *f*.

** Fig. 9.** Evolution of the profiles for a scenario with *f=1, k=1/2, m=1/2*.

a) Evolution of the aperture width. The zone of fast widening has penetrated deeper along the fracture with respect to the isolated case;

b) Evolution of the pressure head distribution. There is a jump depicting the decrease of the pressure head at Node 1;

c) Concentration along the fracture in units of the equilibrium concentration *c eq* . The zone of low concentration is penetrated deeper with the respect to the isolated case;

d) Evolution of the flow rate. The flow through the different fractures is not equal. In the initial moment the flow rate through Fractures 1 and 2 is equal. Afterwards with the widening of Fracture 1, the flow rate through Fracture 1 increases faster. The flow rate through Fracture 2 is increasing relatively fast until, because of the widening of Fracture 0, the pressure head at the Node 1 becomes close to the maximal one at the entrance. From this moment on the flow rate through fracture 2 remains almost constant.

The color code is the same for a), b) and c) and is depicted on Fig. 9c. BT=555 years

We increase the value of *f* to *f=1*. This is a symmetrical structure, because the initial values of the aperture widths are the same for all fractures, and their lengths also. Of course, this is true only at the beginning of the simulation. Later on, because of the dissolution along Fracture 0 and Fracture 1, their resistances decrease, while that of Fracture 2 remains constant. With respect to the cases already discussed the widening of Fracture 0 is more effective (Fig. 9a). The opening zone progresses faster towards the exit. The differences in the pressure head distribution are even easier to note. There is a clearly visible change at Node 1 (see Fig. 9b). Due to the presence of Fracture 2 the gradient increases along Fracture 0. The resistance of Fracture 2 is already sufficiently low and the flow rate through it is comparable to the flow rate through Fracture 1 (Fig. 9d). The profiles of the concentrations are depicted in Fig. 9c. Because of the deeper penetration of the widened zone (see Fig. 9a) and the higher flow rate through Fracture 0, the zone of lower concentration has penetrated deeper along the fracture. The reduction of breakthrough time for this scenario is already more than 300 years.

As we have already seen the influence of Fracture 2 is quite important. Until this moment we have worked with values of *f** £ 1* . In other words, the initial aperture width of Fracture 2 was always smaller or equal to the one of the single fracture. Our next set of experiments will study the behavior of the system at high values of *f*.

** Fig. 10.** Evolution of the profiles for a scenario with *f=4.33, k=1/2, m=1/2*.

a) Evolution of the aperture width. Fracture 0 grows independently of Fracture 1 until the moment of breakthrough at 90 years. Fracture 1 remains unchanged during this period. After the breakthrough event of Fracture 0, the development of Fracture 1 starts under similar conditions – highest possible hydraulic gradient, and low input concentration;

b) Evolution of the pressure head distribution. The pressure head distribution resembles the evolution of two fractures – shifted in time. The second one starts to evolve after the jump in the value of the head at Node 1 – 90 years;

c) Concentration along the fracture with respect to the equilibrium concentration *c eq* . The evolution of the concentration profiles also shows two independently evolving fractures;

d) Evolution of the flow rate. Because of the high initial aperture width of fracture 2, the initial flow through Fracture 0 is also high. The contribution from the flow through Fracture 1 is negligible. After the breakthrough of Fracture 0, the flow through Fracture 2 remains constant. The flow through Fracture 0 remains close to it until the moment of breakthrough of Fracture 1, which marks the breakthrough for the whole system.

The color code is the same for a), b) and c) and is depicted in Fig. 10b. BT=221 years

The arrows indicate restricted breakthrough of Fracture 0, where the flow rates are limited by the finite resistance of Fractures 1 and 2.

For *f=4.33*, the results are depicted in Fig. 10a-d. The influence of Fracture 2 on the pressure head at node 1 is significant (Fig. 10b). The value there is close to the minimal value of the exit boundary condition. This is the reason for the separate evolution of Fractures 0 and 1. There is almost no change in the aperture widths of Fracture 1 during the first 90 years of evolution of Fracture 0 (before breakthrough) (Fig. 10a). The jump at 90 years marks the breakthrough of Fracture 0. This is also visible from the concentration profiles (Fig. 10c), the pressure head distribution (Fig. 10b), and the jump in the flow rates in Fig. 10d. After this, the value of the pressure head at Node 1, changes from 0 m to 100 m – the maximal possible value. At that moment dissolution along Fracture 1 becomes effective. The flow through Fracture 2 remains constant, because there is no dissolution along this conduit. The evolution of Fracture 1 takes roughly the same time, and a second breakthrough event for the system is depicted by the abrupt changes in the profiles (see Fig. 10 a-d) at 221 years.

As a last example we increase *f* to 433.33. This case is surprisingly interesting. Normally one should expect that a further increase of the initial aperture width would result in further decrease of the breakthrough time of the system. This is not the case. The breakthrough time increases by almost 100 years. This “simple” system of three fractures (see Fig. 7a) is actually not so simple. The profiles are depicted on Fig. 11a-d. At a first glance the only change is in distribution of for the pressure head (Fig. 11b). A closer look shows that the aperture width of Fracture 0 after its local breakthrough event is several times larger than the one for the previous simulation. This is easy to explain. The limitations for further evolution of Fracture 0, in the previous scenarios, are Fractures 2 and 1. Now, Fracture 2 is so wide, that it does not limit flow through Fracture 0 and allows considerable opening even before the breakthrough of Fracture 1. After breakthrough of Fracture 0 the head at the junction increases and the evolution to breakthrough of Fracture 1 is initiated. Since the increase in head is smooth now, its evolution is delayed compared to the previous case.

** Fig. 11.** Evolution of the profiles for a scenario with *f=433.33, k=1/2, m=1/2*.

a) Evolution of the aperture width. Fracture 0 grows independently of Fracture 1 until the moment of breakthrough at 90 years. Fracture 1 is unchanged during this period. After the breakthrough event of Fracture 0, the development of Fracture 1 starts under similar conditions – highest possible hydraulic gradient, and low input concentration. The new point with respect to the previous scenario is the further development of Fracture 0 because Fracture 2 is not limiting at this moment;

b) Evolution of the pressure head distribution. The pressure head distribution resembles the evolution of two fractures – shifted in time. The second one starts to evolve after the increase in the value of the head at Node 1. In this case this is a slow increase;

c) Concentration along the fracture in units of the equilibrium concentration *c eq* . The evolution of the concentration profiles shows two independently evolving fractures;

d) Evolution of the flow rate. Because of the high initial aperture width of fracture 2, the initial flow through Fracture 0 is also high. The contribution from flow through Fracture 1 is negligible. After the breakthrough of Fracture 0, the flow through Fracture 2 continues to increase together with the development of Fracture 0. The flow through Fracture 0 remains close to it until the moment of the breakthrough of Fracture 1, which marks the breakthrough for the whole system.

The color code is the same for a), b) and c) and is depicted on Fig. 11c. BT=221 years

Fig. 12 depicts the influence of the parameter *f* on the breakthrough time of the system.

As long the aperture width of Fracture 2 is small, *f*<0.1 little flow is delivered from Fracture 0 and the breakthrough time is close to that of the isolated fracture. With increasing 0.1<*f*<10 the breakthrough time of the system drops. The minimum breakthrough time would occur when two conditions are valid. First the aperture width of Fracture 2 is sufficiently wide such that initially the total head drop is along Fracture 0. Then breakthrough time of Fracture 0 is that of an isolated fracture with length *L/2*. Secondly *f* is sufficiently small to allow a significant head at the junction afer the breakthrough. Then breakthrough of Fracture 1 is also close to that of an isolated fracture with length *L/2*. The total breakthrough time of the system then is the sum of both and an upper limit can be derived from eqn. 8 as Tb (fmin)/Tb single=0.32. This is close to the value at *f*=10 in Fig. 12.

** Fig. 12.** The influence of the parameter f on the breakthrough time of the system from Fig. 7a. – numerical calculations. See text for discussion.

When *f* increases further the breakthrough time increases. For *f* ® ¥ all flow from Fracture 0 leaves the system through Fracture 2. Therefore Fracture 1 receives no flow and widening along it is excluded, and the breakthrough time of the system becomes infinite.

** 3.2. The influence of the position of the exchange node**

Fig. 13 depicts the influence of the parameter *k* on the breakthrough time of the system. We select several values for *f*, representing interesting regions of the breakthrough curve from Fig. 12. These are: *f=0.33, f=1, f=4.33, f=433.33*. Changing the value of the parameter* k* corresponds to a movement of the position of Node 1 along the fracture. It is clear that the maximum of the breakthrough time will be at values of *k = 0* or *1*. In these cases the system develops like a single conduit. With the movement of the junction along the fracture, the breakthrough times decrease. The position of the minimum depends on the value of *f* and *m*. For large values of *f* and *m=1/2* the minimum is *k=1/2*. Then both fractures develop in sequence and independently of each other. The breakthrough time of the system depends only of the lengths of Fractures 0 and 1. The minimum is attained for *L 0=L 1*, i. e. *k=1/2*.

** Fig. 13.** Influence of the position of the junction to the evolution of the system (*m*=1/2).

The minimum of the curves is changing with the change of the *f*. For values of *f* much greater than 1, the minimum in the breakthrough time is when the junction is in the middle of the single fracture. With decreasing *f*, the minimum moves towards the entrance of the fracture.

In the following we will discuss the case with *k=m=1/2*. Depending on the value of *f* two regions of the impact of exchange flow can be envisaged. For small *f<<1* exchange flow is small with respect to the initial flow through Fractures 0 and 1, because of the resistance of Fracture 2. Note that this resistance increases with the third power of *f*. In the extreme case *f=0* the system behaves like an isolated fracture. When the resistances of all three fractures are of similar magnitude, i.e. *f=1*, then the existence of exchange flow acts as positive perturbation to reduce breakthrough times. This is the situation we have encountered in the system of Fig. 8. (*f=0.33*) and Fig. 9 (*f=1*), where slight reductions of breakthrough times have been observed.

Now we turn to the case *f>1*. If the resistance of Fracture 2 is much smaller than that of Fracture 0, the latter fracture experiences a high drop of the hydraulic head. Under these conditions it behaves like an isolated fracture with a head drop close to *H* and experiences breakthrough accordingly. After breakthrough the head at node 1 becomes close to *H*. During this time Fracture 1 has increased its aperture width slightly by according to the small head, which has acted along it. After breakthrough it experiences the full hydraulic head and breakthrough occurs like for an isolated fracture with aperture width ( ). The total breakthrough time for the system is just the sum of these of the two subsequent events. Such scenarios are shown by Fig. 10 (*f=4.33*) and Fig. 11 (*f=433.33*). In both cases the total breakthrough times are equal since *R 2* is larger by a factor of 81 or a factor of 81 × 10 6, respectively. In other words, in both cases *R 2<<R 0, R 1*.

This behavior is also reflected by Fig. 12, which shows the dependence of the total breakthrough time on *f*. For *f>2* we find an almost constant but significant reduction of total breakthrough time. For *f<0.33* breakthrough times are those of the isolated series of Fracture 0 and Fracture 1. The intermediate range of perturbation is between *0.33<f<2*.

## 4. Fracture embedded into a two-dimensional network

We have already shown the importance of the flow lost from the fracture. The modeling setup was as simple as possible – only one point of outflow along the length of the conduit. The next step is to embed the conduit into a two-dimensional fine network. By this way flow can be exchanged at many nodes.

** Fig. 14.** Geological settings modelled by our computer simulations.

Fig. 14 shows the geological setting used as a model for the computer simulations. It depicts a limestone terrain, with a system of small fractures. Embedded into this is a central fracture with larger initial aperture width. The river on the top of the limestone bed supports a constant head boundary condition at *H [m]*. The right hand side of the block is open to flow at constant head at *h=0*.

To model this geological setting, we create the following idealized structure – Fig. 15. The fine fissure system is represented by a two dimensional 100x51 network of fractures. The fractures parallel to the central conduit are termed horizontal fractures. The others – vertical. A constant head boundary condition with head *H* at the left hand side input boundary and with *h=0* at the right hand side output is typical for the early evolution of karst. The upper and lower boundaries are impervious. Only a small amount of flow is carried through these fractures. Dissolution is active along them and they are widened continuously. This lowers their resistance, and the flow increases. After breakthrough the flow continues to increase, until there is not enough water on the input site to support the constant head boundary condition. At this moment the computer runs for this study were terminated. They could be continued further, but then a constant recharge boundary condition must be used.

** Fig. 15.** Modelling domain for the geological settings depicted by Fig. 14. The red line depicts the central wide conduit (initial aperture width (*A 0*). The initial aperture widths of the fine fractures are *a 0* (see Table 2).

** Table 2 **

## 4.1. Evolution of the central channel (Scenario A)

We start with the first experiment. The head on the left hand side is fixed at 100 m. The aperture width of the central channel *A 0* is 0.03 cm, and the aperture widths of the fine net fissures are *a 0*=0.02 cm (uniform net). The results are depicted in Fig. 16. It shows the following profiles along the central fracture:

- Pressure head distribution – Fig. 16a;
- Flow through the central channel – Fig. 16b;
- Aperture widths along the channel – Fig. 16c;
- Concentration along the central conduit, depicted as the ratio between the actual concentration and the equilibrium concentration – Fig. 16d.

We compare these results with those for the isolated fracture from Fig. 5. They are depicted in Fig. 17. The colors of the curves on both figures show the regime of the flow through the conduits. Red lines depict the evolution under laminar flow conditions. Green is used to outline the transition between laminar and turbulent flow, and actually depicts the situation shortly before and shortly after breakthrough. The evolution under the turbulent flow regime is depicted by the blue curves. The breakthrough time for the isolated conduit is about 800 years, and only 86 years for our current scenario.

** Fig. 16.** Evolution of the profiles along a non-isolated conduit for: a) pressure head (the numbers denote the sequence in time); b) flow rate; c) aperture width; d) concentration relative to *c eq. *Red lines – 0 to 80 years every 10 years;

Green lines – 82 to 90 years every 2 years;

Blue lines – 90 to 185 years every 20 years.

** Fig. 17.** Evolution of the profiles along an isolated conduit for: a) pressure head (the numbers denote the sequence in time); b) flow rate; c) aperture width; d) concentration relative to *c eq. *Red lines – 0 to 700 years at every 100 years;

Green lines – 710 to 740 years at every 10 years;

Blue lines – 750 to 850 years at every 10 years

There is one very important point. Initially the pressure head is evenly distributed and linearly decreasing along the isolated and the non-isolated channel, and also along the entire fine network (Figs. 16a and 17a). This is the reason to select the uniform initial distribution of the aperture widths in the fractures of the fine network. There is no head difference perpendicular to the impervious boundaries. Consequently there is no exchange of flow between the central conduit and the fine network. This means that the evolution of both fractures starts from exactly the same point, and all the differences developing in the later stages are due to the presence of the fine network connected to the conduit. These differences are clearly visible even in the earliest stages of the evolution of both channels.

During the first 80 years for the non-isolated case (Fig. 16c) and the first 700 years for the isolated fracture (Fig. 17c) the dissolutional widening is propagating downhead. For the non-isolated fracture, because of the slower widening of the fractures in the fine network (Eq. 6), the pressure difference in the direction to the network is high and the amount of flow leaving the central fracture is also high. This maximum is at the end of the widened zone. From there on part the loss of water into the continuum decreases and reaches zero at the exit of the fracture (Fig. 16b). In contrast, the profiles of the pressure head for the isolated case (see Fig. 17a) are smoother. The curves of the flow rate are horizontal lines because no flow can be lost. The most narrow part of the fracture is acting like a bottleneck and is limiting the flow along the conduit (see Fig. 17b). As it opens flow increases in time.

The profiles of the concentration show a fast exponential increase to a value close to equilibrium, for the isolated case (Fig. 17d). The concentration remains almost constant to the exit. For the non-isolated case the profiles are not so steep. The concentration increases almost linearly along the fracture (Fig. 16d). At the beginning of the unwidened part it reaches a stable value, which remains practically unchanged along the channel, and is lower than that for the single case. The reason for this difference is the higher inflow of aggressive solution into the entrance of the channel.

The transition period and the breakthrough event (the green profiles, 80-90 years for the non-isolated, and 710-740 years for the isolated conduit) start when the widened part reaches the exit (see Figs. 16c and 17c). The pressure head has still the same shape but now the zone of the high head is almost at the end of the fracture for the non-isolated case (Fig. 16a). The profiles for the isolated case show that the zone of the high head has moved towards the exit, but again as in the early stages of the evolution, the profiles are not as steep as in the non-isolated case (Fig. 17a). The flow distribution shows that the zone of the maximal outflow from the fracture to the surrounding environment is now close to the exit of the fracture, following the change in the pressure head distribution (Fig. 16b). Again the concentration profiles along the isolated channel are steeper than the ones along the non-isolated conduit - Figs. 16d and 17d. The dissolution rates at the end of the isolated channel are considerably lower than those at the end of the non-isolated one. This is the reason for much faster widening of the exit zone, and earlier breakthrough for the non-isolated case. The moment of breakthrough is clearly visible on all pictures (86 years for the non-isolated and 740 years for the isolated scenario). After the breakthrough the pressure head distribution changes and is more evenly distributed along the fractures. The zone of the high head has moved backwards along the conduits (Figs. 16a and 17a). At the same time the zone of considerable widening has reached the exit. A dramatic increase (by almost two orders of magnitude) of the flow rates is depicted in Fig. 16b and Fig. 17b. The increase of the concentration now is almost linear along the conduits, for both cases, but the values, at their exits are far away from equilibrium, much lower than the switch concentration (Fig. 16d and Fig. 17d). After the breakthrough, the widening along the fracture stays even along its entire length.

There is an interesting detail on the curve depicting the flow rate along the fracture in the non-isolated scenario at 88 and 90 years. Both curves are not smooth at the end. This is visible close to the exit for the curve depicting the situation at 88 years, and clearly visible further upstream for the curve depicting the flow rate at 90 years. The reason is the growth of an exit fan, into the fine network. This will be discussed later.

After the dramatic increase of the flow rates along the conduits, the flow regime changes from laminar to turbulent. This part of the evolution is depicted by the blue profiles on Fig. 16 and Fig. 17 (100-197 years, for the non – isolated case, and 750 – 850 years, for the single conduit). The zone of the high head continues to move backwards, for the isolated fracture, and the distribution is becoming smoother (Fig. 17a). Because of the even widening with high linear dissolution rates along the entire profile the aperture widths profiles are approaching equal widths and consequently the head distribution becomes linear as at the beginning. The non-isolated conduit shows more complex behavior because of the entrance and the exit fans developing in the fine fractures network. The numbers on the head profiles (green and blue) in Fig. 16 and Fig. 17 indicate their sequence in time. For the isolated conduit, the flow rate increases continuously following the widening of the narrowest part of the fracture (Fig. 17b,c). For the non-isolated case, because of the strong influence of the fine fractures network, the picture is different. There are regions where flow is coming from the net and regions where it is lost to it.

As expected the concentration along the fractures, at this stage of the evolution, is close to zero and the widening is even along them. This is depicted by Figs. 16c, d and Figs. 17c, d.

4.2. Evolution of the 2D network

After having discussed the evolution of the central channel we now turn to the evolution of its environment, the 2D-net. This is discussed for Scenario A where all fractures have equal aperture widths, and for Scenario B, where the net consist of fractures with statistically distributed aperture widths (lognormal, ). Whereas this difference has no significant influence to the evolution of the central channel it shows new features in the net.

** Fig. 18.** Evolution of fracture aperture widths for the standard scenario A. The color code depicts the fracture aperture widths in centimetres. The head distribution is illustrated by the black lines of constant head in steps of 5 m. The numbers in the lower right corner depict the time.

** Fig. 19.** Evolution of flow rates for standard scenario A. The color code depicts the ratio between the flow rate through the current fracture and the maximal flow rate, which occurs at that time in the net. The value of the maximal flow rate in cm 3/s is depicted at the upper right corner. The head distribution is illustrated by the orange lines of constant head in steps of 5m. The numbers in the lower right corner depict the time.

Note that all fractures carrying flow 10000 times smaller than the maximal are not displayed.

The evolution of the aperture widths for both scenarios A and B is presented in Fig. 18 and Fig. 20. The color code depicts the fracture widths in centimeters, and the black lines show the head distribution inside the domain. We expect a similar evolution of both during the period before breakthrough. This is clearly visible in Figs. 18a and Fig. 20a. In both cases a conduit propagates downstream along the central channel. The outflow into the surrounding environment can be visualized from the pressure head distribution. In the statistical case some fissures are initially wider than 0.02 cm. They open relatively fast if the conditions at their entrances and exits are suitable (high pressure gradient, and low calcium concentration for example). This is seen for the widened fractures on the left hand side boundary and for the ones around the tip of the central conduit (see Fig. 20). The flow distribution for Scenario A is illustrated by Fig. 19. The color code depicts the ratio between the flow rate through the actual fracture and the maximal flow rate, which occurs in some fracture of the network. The value of this maximal flow rate is depicted at the upper right corner of every figure. All fractures carrying flow 10000 times smaller than the maximal one are omitted from the figure. The orange lines depict the pressure head distribution as in Fig. 18. Fig. 18a depicts the situation at 70 years. The flow out from the central fracture to the surrounding environment starts close to the left hand side boundary. It increases with the distance from the entrance and becomes maximal at the tip of the opened central conduit. After that it decreases. Close to the exit there is no exchange between the central conduit and the surrounding network. After the tip in the direction of the right hand side border, there is a zone, where the pressure lines are perpendicular to the central fracture. This means that the flow, which leaves the tip is diffused into a large area of fractures. Consequently the aggressive water widens the fractures closest to the central one. The more distant fractures remain largely unaffected because of the large area and the relatively high calcium concentration. The flow distribution on the right hand side boundary is even and still unaffected by the penetrating channel.

** Fig. 20.** Evolution of fracture aperture widths for the standard scenario B. The color code depicts the fracture aperture widths in centimetres. The head distribution is illustrated by the black lines of constant head in steps of 5 m. The numbers in the lower right corner depict the time.

Fig. 18b and Fig. 20b illustrate the situation after 85 years for the uniform case and after 81 years for the statistical one. This is shortly before breakthrough. The significantly widened part of the central conduit is already close to the exit. The main outflow is at the tip of the widened part. This is also visible at Fig. 19b. There is an important difference compared to the situation at 70 years. The zone of the outflow is close to the right hand side boundary. The main consequence is that all the fluid from of the central fracture is not diffused into a large area of unwidened fractures. It is directly connected to the open flow boundary and the water is channelled through the neighbouring fractures parallel to the central conduit. This is the shortest path out of the domain. The central conduit carries the main outflow at the right hand boundary. The outflow decreases in the neighbouring fractures and after some distance it becomes constant. The result is visible on Fig. 18b. Some fractures starting from the central conduit and perpendicular to it are wider than the rest of the fine network. In the statistical case (Fig. 20b) the number of the widened fractures is higher and their distribution is not symmetric. But the main result is similar: Fractures start to develop between the tip of the widened region of the central conduit and the right hand side open flow boundary condition. They belong to the shortest pathways to the exit.

Figs. 18c and 20c depict the situation shortly after breakthrough – for both scenarios. The flow through the central conduit is already turbulent. The narrow part at the end of the central fracture is removed and the head distribution is changed. For the uniform net (Fig. 18c) the pressure lines move backward the gradients are more evenly distributed after 110 years, Fig. 18d. There is a fan at the downstream end of the central fracture. The flow through it is laminar, because the exits of the channels still have relatively small aperture widths. This is the reason for the head difference between the fractures belonging to the fan and the central conduit. A small amount of water between them flows mainly in the direction of the central conduit (Fig. 18c, Fig. 19c, Fig 20b). The exit fractures of the fan are widened until a breakthrough event happens there, and the flow through the whole fan becomes turbulent. This is what is observed in both cases (Fig 18d, Fig. 20c). The flow through the fan becomes turbulent and the pressure difference between its fractures and the central conduit is close to zero. The main inflow to the fan region has moved backwards. Some perpendicular fractures start to develop on both sides of the central conduit. The water entering them is still sufficiently aggressive for relatively fast widening. This increases inflow from the central conduit. The water is not diffused into the rest of the network, because of the exit fan. The influence of the right hand side boundary condition propagates upstream deeper inside the domain and the same mechanism, fast widening of the fractures comprising the shortest pathway to the developed part of the fan, is repeated. Consequently the fan is propagating upstream, on both sides of the central fracture.

This is observed in Figs 18d and 20c. The flow through the fan is already turbulent as depicted by the red, thick fractures. There are two regions of the head distribution. The first region starts at the left hand side of the domain and continues to the region of the fan. The flow is mainly directed from the central region into the unwidened fractures surrounding it as can be seen from the pressure lines at their tips. This means that the outflow is from the central fracture and also from the closest ones parallel to it. The water entering the unwidened part of the domain has lost most of its dissolutional power in the area where the fan grows actively. Therefore it cannot significantly alter the net outside.

The second region of the pressure distribution starts at the upstream end of the exit fan and propagates to the end of the network. Here the situation is exactly the opposite. The pressure lines are bended backwards, showing inflow from the outer unwidened part into the fan and into the central fracture. The concentration of the water coming from the fine fractures is close to saturation and cannot influence the evolution of the widened part. Consequently, the widening of the fan region is limited, and the only direction of growing is upstream towards the entrance of the domain. This is observed at the later stages of the evolution, as shown for the statistical case in Fig. 20d. The exit fan is wider and not symmetrical, but the same behavior is observed. It stops to grow in the direction perpendicular to the central conduit. Due to more favorable pathways, some channels grow from the entrance (left hand side boundary) downstream towards the exit fan.

The situation 30 years later is depicted in Fig. 18e for the uniform and in Fig. 20e for the statistical scenario. In the uniform net the fan is growing further upstream towards the entrance. As mentioned above there is no growth in the direction of the impervious boundaries, because of the saturated inflow coming from the fine network. The outflow is concentrated to the fan and the central fracture, (see Fig. 19e). In contrast to the right hand side, the changes in the left hand region are significant. There is a second fan, propagating down stream. It is formed by the central fracture and the two parallel neighboring fractures. In contrast to the exit fan it is growing from the neighboring fractures in the direction of the central conduit (see Fig. 18e). The growth of the exit fan is supported by the relatively aggressive solution flowing out of the central fracture into the fine network.

The situation for the fan at the entrance is different. The horizontal fractures at the left hand side border of the domain are widened. The initial aperture width (*a 0*) of these fractures is 0.02 cm for this scenario. The breakthrough time of the single fracture is proportional to the 3 rd power of (1/a 0), for n=4 (see Equation 8). A fracture with 0.02 cm initial aperture width, will take almost 3.4 times longer time for breakthrough to the exit than a fracture with initial aperture width 0.03 cm. A period of 140 years is long enough to observe widening in the fractures forming the entry part of the fine network. During the evolution to breakthrough, the fractures distant from the central conduit (in the direction of the impervious boundaries) have experienced different pressure head at their exits than the fractures close to it during the whole evolution so far (see Figs. 18a-d). This explains the smaller aperture widths of the entry fractures close to the central conduit.

The entry fan starts to grow from the first pair of fractures neighboring the central conduit and parallel to it (Fig. 18e). All entry fractures (including the first pair) have enough time to widen considerably. This increases the pressure head at the nodes of their exits. Consequently the pressure distribution along the entire entrance part of the domain is changed. There is a zone of inflow to the central conduit. It starts at the ends of the fractures forming the left boundary of the domain, and ends close to the tip of the entry fan. Part of the water is diverted directly to the central conduit, and the other part continues to flow parallel to it. This water is still aggressive, and causes fast widening. The presence of the exit fan diverts a relatively high amount of water in the direction parallel to the central conduit, and consequently the growth of the entry fan is fast in both directions (forward – parallel to the flow), and in the direction perpendicular to the central channel. Actually there is no mechanism to stop its growth. The water, entering the domain on the left hand side is aggressive, and if there is a wider central conduit, then sooner or later the evolution will come to a phase when the flow will be diverted towards this wider fracture. At this moment the entry fan starts to grow until the whole domain is conquered. The pressure distribution, in the transition zone between the entry and the exit fans, matches the switch between the mechanisms governing their growth. Next to the end of the entrance fan, the flow is already parallel along the entire central fracture, and the exchange is close to zero. This depicts the beginning of the zone, where the exit fan is developing.

Fig. 20e depicts the situation for the statistical scenario at 140 years. The widened fractures on the left hand side are forming the entrance fan, which is already connected with the exit fan. The whole central zone exhibiting widened aperture widths (red region) can be considered as equivalent to one huge central fracture. As already discussed, there is no mechanism limiting the growth of the entry fan. It is extending further in horizontal and in vertical direction at both sides of the central zone. The area of the inflow to the central region is extending from the left hand side of the network to the right border . There is no longer flow out from the central channel to the surrounding network. Consequently development of the exit fan is stopped.

In contrast to the uniform case, only some of the horizontal entrance fractures are widened (see Fig. 20e). The others remain unchanged. On the other hand, the ones, which are widened extend deeper. The reason is the statistical distribution of the initial aperture widths of the fractures. Any of the widened and deeper penetrating channels, can be considered as a central wider fracture for the region, where it grows. All the conclusions, so far, are true also for such a small sub domain. The widened path continues to grow in the direction of the highest hydraulic gradient. These pathways can reach directly the right hand side boundary or can divert their growth in the direction of the widened area in the central part of the domain, and connect to it. Entrance and exit fans develop on both sides of these pathways. An example for such a place is depicted by the black rectangle on the figure. Actually, when these pathways connect either to the central widened zone or to the right hand side boundary they exhibit their local breakthrough event.

The situation after 180 years of evolution is depicted in Fig. 18f for the uniform scenario and in Fig. 20f for the statistical one. As expected, the entry and the exit fans are already connected in the uniform scenario. The entry fan continues to grow. The pressure distribution is similar to the one for the statistical case at 140 years. The flow is directed towards the huge widened zone around the central conduit (Fig. 19f) and there is no area of outflow from the central fracture to the fine network. As already discussed, the consequence of this is the end of the evolution of the exit fan. We expect continuous growth of the entry fan in both directions – horizontal, to the exit of the domain, and vertical, in the direction of the impervious borders. This is depicted in Fig. 20f – the statistical scenario, where almost the whole domain is conquered by the entrance fan.

## 4.3. Influence of the exchange flow on the breakthrough time

Our next goal is to investigate the influence of the hydraulic conductivities of the fine fracture network on the evolution of the central fracture.

** Fig. 21.** Breakthrough time for Scenario A (squares) and Scenario B (triangles), as a function of initial aperture widths of the continuum fractures *a 0*.

To this end the initial aperture widths *a 0* of the fine fractures are varied from 10 -5 to 0.03 cm and *A 0* is fixed at 0.03 cm for the central conduit (see Table 2). Everything else is left unchanged. The dependence of breakthrough times on the initial aperture width *a 0* for both scenarios A and B is shown on Fig. 21. The two curves look almost identical, which means that for the period until breakthrough, the evolution of the two systems (uniform and statistical), is similar. The breakthrough time is maximal when:

- The aperture widths of the fractures of the surrounding network are negligibly small (
*A 0>>a 0*); - The aperture widths of the fine fractures are the same as the aperture width of the central fracture (
*A 0=a 0*). This is possible only for the uniform scenario A.

In the first case the resistances of the surrounding fractures are so high, that the flow they can carry is close to zero. Consequently the exchange between the central fracture and the surrounding network is also close to zero and the central fracture can be regarded as an isolated one.

In the other extreme case, the head distribution is uniform along the whole network (only for the uniform networks). There is no pressure difference between the central fracture and the surrounding network. All flow lines are parallel to the upper and to the lower boundaries and there is no exchange of water. In both cases, there is no flow out from the central fracture into the continuum and the breakthrough times should be equal to that of an isolated fracture. This is exactly what can be seen on the picture.

For finite aperture widths a 0<0.03 cm flow from the central fracture into the net exists. Consequently by the higher inflow of aggressive water into this fracture the dissolutional widening at the exit is enhanced and breakthrough times become shorter, with increasing *a 0*.

A reasonable question is: What would happen if the initial widths of the surrounding fractures were increased even further, beyond the size of the central one. Then their breakthrough times would be shorter in comparison to the breakthrough time of the central fracture, and the breakthrough event will happen in the surrounding network. Our goal, however, is to show the evolution of the central conduit. Therefore we will not continue the curve in Fig. 21 further.

* 4.4. Extended scenarios *

In all cases discussed so far, the pressure head on the input site of the domain is extremely high (*H=100 m*). At the same time the initial aperture widths of the fine fractures have values relatively close to the ones of the central conduit:*a 0=0.02 cm* and *A 0=0.03 cm*, respectively.

** Fig. 22.** Evolution of the fracture aperture widths for extended scenario A. Everything is the same as in the standard scenario A (see Table 2.1.1), except the hydraulic head at the input side. In this case it is *H=10 m*. The color code depicts the fracture aperture widths in centimetres. The number at the lower right corner depicts the time.

We start the next set of simulations with a change of the pressure head applied to the input site. It is fixed at *H=10 meters*. The head at the output site remains unchanged – *h=0 meters*. The other initial and boundary conditions are the same as for the standard scenarios A and B (see Table 2). Fig. 22a-d depicts the evolution of the fracture aperture widths for the uniform scenario. Lowering the hydraulic head by one order of magnitude increases the breakthrough time of a single conduit (see Equation 8). The breakthrough time of the central conduit in the case of *a 0*=10 -5 cm when exchange flow is negligible is about 800 years under the pressure head of *H=100 meters*. With a 10 times lower head it has to be in the range of 17000 - 18000 years. The breakthrough time of our uniform scenario under the pressure head *H=10 meters* is close to 7500 years. Even under the lower hydraulic gradients, the breakthrough time of the aquifer is reduced due to the presence of exchange flow. Fig. 22a shows that the central fracture starts to develop downstream. The exit fan starts to grow and at 7800 years (shortly after the breakthrough of the central conduit) some of the fractures are considerably wide – Fig. 22b. Fig. 22c depicts the situation after 9100 years. The exit fan has migrated backwards. At the same time an entry fan starts to grow downstream towards the right hand side boundary. After 9400 years, the evolution of the exit fan is stopped, but the entry fan continues to penetrate deeper into the block (Fig. 22d). After this moment the constant head boundary condition cannot be supported anymore and the simulation is terminated. The difference to the standard scenario is in the penetration distances of both fans. These are much shorter. This is an important observation, because it shows that the difference between the time scales for the evolution of both hydraulic systems (the fine network and the central conduit) is important. The fine network evolves slower with respect to the central fracture, than for the standard scenario A with *H*=100 m.

Fig. 23 a-d depicts the evolution of the statistical case. The results are similar to those for the uniform scenario.

** Fig. 23.** Evolution of the fracture aperture widths for extended scenario B

Everything is the same as in the standard scenario B (see Table 2), except the hydraulic head at the input side. In this case it is *H=10 m*. The color code depicts the fracture aperture widths in centimetres. The number at the lower right corner depicts the time.

In our next experiment we reduce the head at the input to 3 meters. This value is realistic for natural karst aquifers. Everything else is like in the standard scenario A (see Table 2). Figs. 24a-d depicts the evolution of the fracture aperture widths for the uniform scenario. The differences in the time scales for the evolution of the fine network and the central fracture are already so large, that there is no time for considerable changes in the network. After the breakthrough event, the central conduit continues to evolve under turbulent flow conditions. The concentration along its entire length is low and the rates of widening are close to the maximal ones.

** Fig. 24** Evolution of the fracture aperture widths for extended scenario A Everything is the same as in the standard scenario A (see Table 2), except the hydraulic head at the input side. In this case it is *H=3 m*. The bar code depicts the aperture widths of the fractures in cm. All fractures, which are smaller then 0.08 cm, are omitted from the figure The number at the lower right corner depicts the time.

The next experiment has the following initial and boundary conditions. Everything (including the hydraulic head on the entrance) is as for the standard scenario B (see Table 2). The only change is the average initial aperture width of the fractures in the fine network. In this case it is set to cm. The result is depicted by Fig. 25 a-d. Because of the smaller initial aperture widths compared to scenario B with *a 0*=0.03 cm, the time for the evolution of the fine network is increased. Similar to the low hydraulic gradient case (*H=3m*), there is no time for considerable changes in the fine fractures. This is a confirmation for the conclusion, that when the time scales of both hydraulic systems are similar, then the changes in the fine network are important for the development of the karst aquifer. On the other hand, if these time scales are sufficiently different, then both systems develop almost independently.

** Fig. 25.** Evolution of the fracture aperture widths for extended scenario B. Everything is the same as in the standard scenario B (see Table 2), except the initial aperture widths of the fine fractures. In this case they are statistical distributed with average at *a 0=0.01 cm. *The bar code depicts the aperture widths of the fractures in cm. The number at the lower right corner depicts the time.

The following statements can summarize these few examples of aquifer evolution. The evolution proceeds in two distinct steps. In the first period until breakthrough of the central fracture, flow always is injected from the wider fracture into the net. The main evolution is therefore widening of the fracture, whereas the fissures in the net are not affected. After breakthrough the central fracture widens evenly and continuously about 0.1 cm/year and the hydraulic gradient becomes constant along it. At the same time exit and entrance fans are generated. The first step of entrance fan evolution is a breakthrough event from the input point of the nearest neighbour longitudinal fissure to the next node and then down the vertical fissure connecting to the central conduit. The head difference along this path is *h/j*, where *j* is the number of nodes in the central fracture. For a significant evolution of a fan system about *j/2* of such subsequent events are necessary. By use of eqn. 8 and multiplying the thus obtained breakthrough time by *j/2* one obtains a crude estimation of fan evolution time *T fan* as:

where *T c* is the breakthrough time of the central channel, *j=100*.

If *T fan* is larger than the time elapsed after the breakthrough until the constant head can no longer be supported fans cannot evolve. This is the case for long breakthrough times as in natural karstification with small fracture aperture widths and low hydraulic gradients. For high hydraulic heads, as supported close to dam sites, breakthrough times are short and extended fans can evolve below the dam site as shown recently by Romanov et al. (2003).

## Conclusion

We have discussed the changes of conduit evolution of a single isolated fracture, under constant head boundary conditions with highly undersaturated input solutions, when one allows exchange of flow, either into a tube connected to it, or by embedding the single fracture into a net of narrow fissures with smaller aperture widths, into which water can be transferred. Essentially both means to add a second dimension to the one-dimensional single tube. In all cases when exchange flow leaves the conduit more aggressive solution enters into its entrance and penetration length of dissolutional activity increases. By this way dissolution rates at the exit are enhanced and the time needed until breakthrough becomes shorter. When the fracture is embedded into a network of narrow fissures it looses flow into the net during the entire period of initial evolution until breakthrough. After breakthrough the concentration along the entire conduit becomes low and these highly aggressive waters are lost to neighboring fractures, causing small breakthrough events through them. This creates an exit fan, which propagates upstream. With some delay an entrance fan develops from the input propagating downstream and invading the entire net to a finally maze like structure. This is a result of highly aggressive waters flowing from fractures close to the central channel into it, thereby creating breakthrough. The structure extends all over the net, because once a neighboring fracture has become sufficiently wide it acts like the central one to the outer fissure next to it and so on. The evolution proceeds in two steps. In the first until breakthrough of the central channel, flow is directed into the net and dissolution occurs mainly in the channel. After breakthrough the channel widens evenly with a rate of about 1mm/year. At the same time an entrance fan is created by subsequent breakthrough events from neighboring fissures to the channel. If the time for the evolution of this fan is shorter than the time, when a constant head can no longer be supported, extended maze like structures evolve. In this case correct modeling must consider dissolutional widening in the net. In cases where fans cannot develop the evolution of the aquifer affects mainly the central channels and models neglecting dissolution in the net (Bauer et al., 2003) are a satisfactory approximation.

## Bibliography

Bauer S. 2002. Simulation of the genesis of karst aquifers in carbonate rocks, Tuebinger geowossenschaftliche arbeiten (TGA), C62, 2002.

Bauer S., Liedl, R., Sauter, M. 2003. Modeling of karst aquifer genesis: Influence of exchange flow. Water Resources Research., 39 (10), 1285, 2003.

Beek W.J. and Muttzall K.M.K. 1975. Transport Phenomena. New York: Wiley.

Bear and Verruijt. 1987. Modeling groundwater flow and pollution. Dordrecht, Holland: D. Reidel Publishing company.

Clemens T., Hückinghaus D., Sauter M., Liedl R., Teutsch G. 1996. A combined continuum and discrete network reactive transport model for simulation of karst development. In Calibration and realiability in groundwater modelling – proceedings of the ModelCARE 96 conference held at Golden, Colorado, September 1996, IAHS Publ. no. 237, 309-318 .

Clemens T., Hückinghaus D., Sauter M., Liedl R., Teutsch G. 1997a. Modelling the genesis of karst aquifer systems using a coupled reactive network model. In: Hard Rock Geosciences. IAHS Publ., 241, Colorado: 3-10.

Clemens T., Hückinghaus D., Sauter M., Liedl R., Teutsch G. 1997b. Simulation of the evolution of maze caves. In proceedings of the 12 th International congress of speleology, La Chaux-de-Fonds, Switzerland, Volume 2, pp65-68.

Dreybrodt W. 1988. Processes in karst systems - Physics, Chemistry and Geology. Springer Series in Physical Environments 4, Berlin- New York: Springer, 188p.

Dreybrodt W. 1990. The role of dissolution kinetics in the development of karstification in limestone: A model simulation of karst evolution. Journal of Geology 98, 639-655.

Dreybrodt W. 1996. Principles of early development of karst conduits under natural and man-made conditions revealed by mathematical analysis of numerical models: Water Resources Research 32, 2923-2935.

Dreybrodt W. and Siemers J. 2000. Cave evolution on two-dimensional networks of primary fractures in limestone: In: Klimchouk, A., Ford, D.C., Palmer, A.N., and Dreybrodt, W., (Eds.), Speleogenesis: Evolution of karst aquifers. Nat. Speleol. Soc., USA, 201-211.

Dreybrodt W., Romanov D. and Gabrovšek F. 2001. Karstification below dam sites: a model of increasing leakage. In: Beck, B.F., Herring, J. G., (Eds.). Geotechnical and environmental applications of karst geology and hydrology: Proceedings of the 8 th Multidisciplinary Conference on Sinkholes and the Engineering and Environmental Impacts of Karsts, Louisville, Kentucky, Balkema, p. 131-137.

Eisenlohr L., Meteva K., Gabrovšek F., Dreybrodt W. 1999. The inhibiting action of intrinsic impurities in natural calcium carbonate minerals to their dissolution kinetics in aqueous H 2O-CO 2 solutions. Geochimica et Cosmochimica Acta 63, 989-1002.

Ford D.C. Williams P.W. 1989. Karst geomorphology and hydrology. Unwin Hyman, London

Gabrovšek F. 2000. Evolution of early karst aquifers: From simple principles to complex models: Zalozba ZRC, Ljubljana, Slovenia, 150p.

Gabrovsek F. and Dreybrodt W. 2000a. The role of mixing corrosion in calcite aggressive H 2O-CO 2-CaCO 3 solutions in the early evolution of karst aquifers. Water Resources Research 36: 1179-1188.

Gabrovsek F. and Dreybrodt W. 2000b. A model of early evolution of karst conduits affected by subterranean *CO 2* sources. Environmental Geology 39:531-543, 2000.

Gabrovsek F. and Dreybrodt W. 2001. A model of the early evolution of karst aquifers in limestone in the dimensions of length and depth. Journal of Hydrology 240: 206-224.

Groves C., Howard A. 1994.Early development of karst systems 1. preferential path enlargement under laminar flow. Water Resources Research 30, 2837-2846,

Howard A.D. and Groves C.G. 1995. Early development of karst systems 2. Turbulent flow. Water Resources Research 31, 19-26.

Kaufmann G. and Braun J. 1999. Karst aquifer evolution in fractured rocks. Water Resources Research 35, 3223-3238.

Kaufmann G. and Braun J. 2000. Karst aquifer evolution in fractured, porous rocks. Water Resources Research 36, 1381-1391.

Kaufmann G. 2003. A model comparison of karst aquifer evolution for different matrix-flow formulations. Journal of Hydrology 283 (1-4), 281-289.,

Lauritzen O. N. 1992. Modeling the evolution of channel networks in carbonate rocks. In J. Hudson (Ed), ISRM Symposium: Eurock’92, pages 57-62. Thomas Telford.

Lee C. and Farmer I. 1993. Fluid flow in Discontinuous Rocks. London: Chapman & Hall.

Lowe D. 2000. Development of Speleogenetic ideas in the 20 th century: The early nodern approach. In: Klimchouk A., Ford D.C., Palmer A.N., W. Dreybrodt (Eds), Speleogenesis: Evolution of karst aquifers. Nat. Speleol. Soc.,USA.

Palmer A.N. 1988. Solutional enlargement of openings in the vicinity of hydraulic structures in karst regions: In: 2 nd Conference on Environmental Problems in Karst Terranes and Their Solutions. Assoc. of Groundwater Scientists and Engineers Proceedings, USA, 3-15.

Palmer A. N. 1991. The origin and morphology of limestone caves. Geol. Soc. of America Bull. 103, 1-21.

Romanov D., Gabrovsek F., Dreybrodt W. 2003. Dam sites in soluble rocks: a model of increasing leakage by dissolutional widening of fractures beneath the dam. Engeneering Geology 70.

Shaw T. 2000. Views on cave formation before 1900. In: Klimchouk A., Ford D.C., Palmer A.N., W. Dreybrodt (Eds), Speleogenesis: Evolution of karst aquifers. Nat. Speleol. Soc., USA.

Siemers J. 1998. Simulation von Karst-Aquiferen. PhD Thesis, Universitaet Bremen.

Siemers J. and Dreybrodt, W. 1998. Early development of karst aquifers on percolation networks of fractures in limestone. Water Resources Research 34, 409-419.

White W. B. 2000. Development of speleogenetic ideas in the 20 th century: The modern period, 1957 to the present. In: Klimchouk A., Ford D.C., Palmer A.N., W. Dreybrodt (Eds), Speleogenesis: Evolution of karst aquifers. Nat. Speleol. Soc.,USA.