Evolution of karst conduit networks in transition from pressurized flow to free-surface flow

Most of the existing models of speleogenesis are limited to situations where flow in all conduits is pressurized. The feedback between the distribution of hydraulic head and growth of new solution conduits determines the geometry of the resulting conduit network. We present a novel modeling approach that allows a transition from pressurized (pipe) flow to a free-surface (open-channel) flow in evolving discrete conduit networks. It calculates flow, solute transport and dissolution enlargement within each time step and steps through time until a stable flow pattern is established. The flow in each time step is calculated by calling the US Environmental Protection Agency Storm Water Management Model (US Environmental Protection Agency, 2014), which efficiently solves the 1-D Saint-Venant equations in a network of conduits. Two basic scenarios are modeled, a low-dip scenario and a high-dip scenario. In the low-dip scenario a slightly inclined plane is populated with a rectangular grid of solution conduits. The recharge is distributed to randomly selected junctions. The results for the pressurized flow regime resemble those of the existing models. When the network becomes vadose, a stable flow pathway develops along a system of conduits that occupy the lowest positions at their inlet junctions. This depends on the initial diameter and inlet position of a conduit, its total incision in a pressurized regime and its alignment relative to the dip of the plane, which plays important role during the vadose entrenchment. In the high-dip scenario a sub-vertical network with recharge on the top and outflow on the side is modeled. It is used to demonstrate the vertical development of karst due to drawdown of the water table, development of invasion vadose caves during vadose flow diversion and to demonstrate the potential importance of deeply penetrating conductive structures.


Introduction 1.Speleogenetic models: A short history, aims and results
Karst aquifers are among the most prolific water reservoirs.Due to their heterogeneity and anisotropy, their efficient exploitation and protection face many challenges.The role of solution conduits in karst aquifers has been a topic of numerous studies.Estimates show that conduits carry about 99 % of flow within karst aquifers and present efficient transport pathways for potential pollutants (Worthington, 1999).However, we have only limited insight into karst aquifers; the position of conduit systems is largely unknown, except for the parts accessible for human exploration or encountered directly by drilling or indirectly by geophysical techniques.
Speleogenesis (e.g. the evolution of conduit networks in karst aquifers) has been one of the main topics in karst studies of the last century (Ford and Williams, 2007).Many conceptual models of speleogenesis have been proposed based on field observations (Audra et al., 2007;Ford and Ewers, 1978;Audra and Palmer, 2013;Palmer, 1991) and inference from basic principles of flow.However, to gain insight into the processes governing speleogenesis, different physical models have been built and followed by numerical models that are based on the physical and chemical principles of flow, dissolution and transport.
The main objectives of speleogenetic modelling are to test the conceptual models, to determine and evaluate the role of different geological, hydrological and geochemical factors and to find mechanisms that govern the evolution of conduit networks in karst aquifers.Few examples of direct field application have been published (Epting et al., 2009).Ewers (1982) applied hardware (physical) models made from plaster of Paris or salt, and discovered several key mechanisms that were later largely confirmed and extended by numerical models.Numerical modelling of single conduit evolution (Dreybrodt, 1990(Dreybrodt, , 1996;;Palmer, 1991;Dreybrodt and Gabrovsek, 2000) revealed a feed-back mechanism between flow and dissolution rates and stressed the importance of higher order dissolution kinetics (Dreybrodt, 1990;Palmer, 1991;Dreybrodt, 1996;White, 1977) for the evolution of extended conduits.Such kinetics has been proven experimentally for limestone and gypsum (Eisenlohr et al., 1999;Jeschke et al., 2001).More elaborated models of 2D fractures with statistical aperture fields (Hanna and Rajaram, 1998;Szymczak and Ladd, 2011) showed that nonlinear kinetics is not necessary for the evolution of extended patterns of solution conduits.
The initial stage of speleogenesis is characterised by slow enlargement of proto-conduits, which is accelerated by positive feedback between flow and dissolution rate under constant head conditions.Dissolutional widening increases the flow rate along an initial fracture.Then as the flow rate increases, fresh aggressive solution penetrates deeper into the fracture and in turn accelerates widening and flow rates.This feedback mechanism leads to breakthrough, when flow and widening rate increase by several orders of magnitude in a very short time (Dreybrodt, 1990;Palmer, 1991;Dreybrodt, 1996;Dreybrodt and Gabrovsek, 2000).At breakthrough the initiation stage of conduit development ends and the enlargement stage starts.The time needed to reach breakthrough is termed breakthrough time.

Evolution of a discrete network under pressurised flow conditions
Individual fractures have been assembled into fracture networks in order to model patterns of evolving conduit systems (Lauritzen et al., 1992;Groves and Howard, 1994;Siemers and Dreybrodt, 1998;Kaufmann and Braun, 2000;Liedl et al., 2003).A typical benchmark setting emerged out of the Ewers's hardware models.It includes a plane populated with initial proto channels (fractures/tubes) with inputs and outputs at different hydraulic heads.These models revealed the competition between different pathways connecting inputs to outputs, as already observed by Ewers (Ewers, 1982;Ford and Williams, 2007) in the physical model.
To review some of these basic mechanisms, a simple scenario is shown in Fig. 1.It consists of a plane with a rectangular grid of fractures.The boundary conditions are shown on Fig. 1a: the sides of the network are marked geographically N, S, E, W. No-flow conditions are applied on the N and S boundaries.Water enters the network at two inputs, In1 and In2 at the W side, initially at constant head H = 5000 cm.The whole E boundary presents output at H = 0 m.Initial aperture widths of fractures are set to 0.02 cm, except for the fractures along W-E line connecting In 1 to the output boundary, denoted as P1 , which has slightly larger initial aperture (0.03 cm) and evolves faster than P2 (Fig. 1a), which is fed directly by In 2. Figure 1 shows aperture widths as line widths and dissolution rates as line colours; the warmer the colour the higher the rate.Equipotential lines are also shown on Figs.1a-f, which show the network at different time stages, denoted in each panel in units of breakthrough time T B .At 0.99 T B (Fig. 1a) the high head from the input has penetrated along the widened fractures of P1 deep into the network, and suppressed both the hydraulic gradient and growth of P2.
Figure 2 shows the profile of hydraulic head along P1 (dashed) and P2 (full line) at different stages, as denoted by arrows.The gradient between the tip of P1 and outputs increases in time until the breakthrough.
After breakthrough (Fig. 1b), P1 is widened with the maximum dissolution rate along its entire length.It becomes increasingly uniform and so does the hydraulic gradient along it (see Fig. 1b and Fig. 2).The gradient builds up between the high head region along still prebreakthrough ("plugged") P2 and post breakthrough ("released") P1, which triggers the growth of conduits connecting P2 to P1. Grey arrows show some principle directions of growth.
Two "post-breakthrough" scenarios are envisaged: 1) In Fig. 1c and d, the constant head is kept at both inputs.New connections between P2 and P1 evolve, while P2 also grows towards the exit.The network expands along the existing pathways by growth of new bypasses (some are shown by grey arrows) until all possible flow paths evolve (not shown).Of course, all catchments have limits and such conditions cannot last for long.
2) In Fig. 1e and f, the recharge at In1 and In2 is limited to Q max = 500 l/s.In this case the constant head conditions break, when inflow at the input reaches Q max .At 1.5 T B (Fig. 1e), the head at the input of P1 is about 1/5 of h max (see also Fig. 2) and the gradient from P2 towards P1 is high, as In2 is still under maximal head.P2 integrates with P1, but further expansion of network is suppressed as the head along the growing existing pathways decreases in time.The interested reader is referred to a detailed modelling study on the influence of limited discharge upon the resulting distribution of conduit sizes by Hubinger and Birk (2011).
To summarise: In pressurised flow conditions, the evolution of the network starts with competition of pathways connecting inputs to outputs and continues with their integration and expansion until head gradients along un-evolved pathways are high enough for pathways to breakthrough.The evolution is controlled by the feedback mechanism between the distribution of hydraulic head and growth of new conduit pathways.This interplay is affected by many parameters which reflect local hydrology, geology and geochemistry.
Many other scenarios of early speleogenesis have been modelled to study factors such as the role of geochemical conditions and mixing corrosion, exchange flow between the matrix and conduit network, and the role of insoluble rocks in the evolution of conduits (Dreybrodt et al., 2005).Numerical models have been also used to assess increased leakage at dam sites or other hydraulic structures where unnaturally high hydraulic gradients cause short breakthrough time (Dreybrodt, 1996;Romanov et al., 2003;Hiller et al., 2011).
In real situations the available recharge cannot sustain pressurised flow within the evolving network, and the conduits undergo a transition from pressurised to free surface flow conditions.Most accessible cave systems have undergone such a transition.
Though most models have only considered pressurised flow, Annable & Sudicky (1998) and Annable (2003) developed an elaborate model of the evolution of a single partially filled conduit embedded in variably saturated fractured media under laminar flow conditions.The extension of such a model to networks with turbulent flow remains a future challenge.
Here we develop a model that goes beyond the dynamics depicted in Fig. 1 by incorporating the transition to, and further evolution in, a free surface flow regime.

Evolution of karst conduit networks in the vertical dimension
The vertical evolution of karst has been under debate for more than a century, starting with classical concepts of Katzer, Grund, Davis, Swinnerton, Rhoads and Sinacori and others (Palmer, 2007).The Four State Model of Ford and Ewers (1978) elegantly combines these concepts and relates cave geometry to the density of permeable fissures.Gabrovšek and Dreybrodt (2001) and Kaufmann (2003) modelled a 2D vertical cross-section of a karst system to explore the evolution of karst aquifers in the dimension of length and depth (sensu Ford & Ewers (1978).They have shown the important role of water table drawdown in speleogenesis.These models considered dissolution in the phreatic part of an aquifer only and partly modelled the formation of drawdown vadose passages (Ford, 1988;Ford and Williams, 2007).Conceptual models have been developed that hypothesize the diversion of vadose water and formation of invasion vadose systems (Ford, 1988;Ford and Williams, 2007;Palmer, 2007;Audra and Palmer, 2013).However, these conceptual models have not been tested by numerical models.
In the following sections we describe how the model is built and present two basic modelling scenarios, each with several representative cases.We focus on the description of new mechanisms of flow pathway selection and discuss the results in view of the existing conceptual models.

The model set up 2.1 The conceptual approach
Figure 3 shows a conceptual framework for the modelling presented in this work.We assume a plane populated with conduits with water-soluble walls, similar to that in Fig. 1.Water enters the conduit network at selected junctions indicated by arrows in Fig. 3.The direct recharge into a junction is limited either by the elevation of the land surface (h max ) or by the maximal available recharge Q max ; if the hydraulic head is lower than h max , all available recharge (Q max ) enters at the junction, otherwise the hydraulic head at the junction is equal to h max and only part of the available recharge enters the system.A similar hardware model was discussed by Ewers (1982) who used the term Multiple-input Multi-rank scenario.
The basic workflow of the model follows the same scheme as in the models cited above (e.g.(Dreybrodt et al., 2005) and includes the following steps: 1. Define the network of conduits and boundary conditions (water inlets and outlets).We also assume that: 1.The flow does not depend on the dissolved load.
2. Time scales for flow, dissolution and transport can be separated from the timescale for widening, i.e. the evolution goes through a set of stationary states within which the widening is constant.

The calculation of flow
We assume that the network has passed the initial (inception) stage of speleogenesis and that turbulent flow has already been established in the network.The reader is referred to work of Dreybrodt et al. (2005) for early evolution in the laminar flow regime.One-dimensional turbulent flow is considered within all conduits.The flow could be either pressurised or free surface.
Flow in partially filled conduits is described by Saint Venant equations (Dingman, 2002), which are based on depth-averaged conservation of mass and momentum.Several numerical techniques are used to solve them (Dingman, 2002).Our model invokes an open source package Storm Water Management Model (abbreviated SWMM from here on), developed primarily for flow and transport simulation in sewage systems by the US Environmental Protection Agency (EPA, 2014).SWMM solves the set of Saint Venant equations to the desired approximation and accuracy using successive approximations with underrelaxation (Rossman, 2009).Its use for the simulation of flow in conduit dominated karst systems has been demonstrated by several authors (Peterson and Wicks, 2006;Gabrovšek and Peric, 2006;Halihan et al., 1998).The pressurised flow is accounted for by introduction of a fictitious Preissmann slot (Fig. 4 ) at the top of a conduit's cross-section (Cunge and Wegner, 1964).In this way we transform a pressurised pipe to an open channel without considerably changing the hydraulic characteristics and enable use of the same set of equations for both flow regimes.Friction losses in conduits are calculated by the Manning equation where S f is the friction slope, V the flow velocity, R the hydraulic radius (i.e. the ratio between cross-sectional area of flow and wetted perimeter), n the Manning roughness coefficient, here taken in the range 0.01 < n < 0.02, k a correction factor depending on the unit system used.For the metric system, k = 1 m 1/3 /s.By introducing k, n remains dimensionless.
SWMM enables easy construction of an arbitrary conduit network and many additional elements, such as reservoirs, catchments etc., which could be implemented into future upgrades of the models presented here.

Dissolution and transport
Dissolution rates in karst environments are determined by the reaction kinetics at the rockwater interface (surface controlled dissolution), by diffusion transport of ionic species between the water-rock boundary and the bulk solution (transport controlled dissolution), and, in the case of carbonates, by the rate of CO 2 hydration (Kaufmann and Dreybrodt, 2007).
Each of these mechanisms can be rate limiting under certain conditions.
In the early evolution of conduit networks, the water in protoconduits (sub millimetres to few millimetres in size) is close to equilibrium with the mineral being dissolved and dissolution is mostly surface controlled, by higher order kinetics in case of limestone and gypsum.In turbulent flow conditions, for cases discussed in this work, dissolution in limestone is dominantly surface controlled by first order kinetics, if the input solution has low saturation ratio.Some issues related to limestone dissolution rates in turbulent flow still remain open; scalloped walls of limestone caves suggest that transport control might play an important role under turbulent flow conditions as well (Covington, 2014).
For these reasons we simplify the dissolution kinetics by assuming a linear rate law at the rock-water boundary: (2) () where α s is the kinetic constant, c eq is the equilibrium concentration of ionic species of the rock forming mineral and c s their actual concentration at the surface of the mineral.Ions are transported from the surface into the bulk through a Diffusion Boundary Layer (DBL) of thickness ε (Dreybrodt and Buhmann, 1991).The transport rate through the DBL is given by: (3) () .
D is a diffusion coefficient, ε the thickness of the diffusion boundary layer and c the concentration in the bulk solution.Equating Eqs. ( 2) and ( 3) gives an equation for c s and an expression for the effective rates: (5) ( ); α t depends on the thickness, ε, of the DBL, which is related to the thickness, h, of the viscous sub-layer by Schmid's number (Schlichting and Gersten, 2000): where ν is kinematic viscosity and Sc the Schmidt number, which represents the relation between the viscous diffusion rate and mass diffusion rate.The thickness of a viscous layer over a flat wall is given by (Incropera and DeWitt, 2002): where   is viscous shear stress at the wall and  is the water density.
Viscous shear stress is related to the friction slope S f (8) , where g is Earth's gravitational acceleration.Taking the Manning relation (Eq.( 1)) for S f and inserting ( 8) into ( 7), gives: Inserting Eq. ( 9) into Eq.( 6) and further into Eq.( 4), we get an expression for ε and for the transport constant α t : (10 Most cases that we present in this work assume that st   , so that t   .Therefore, the dissolution rates are transport controlled.Usually higher flow rates bring with them stronger mixing, lower bulk concentrations and higher dissolution rates.In most situations, the rule of thumb is: the higher the flow, the higher the dissolution rate.
The ions entering the water increase its saturation state with respect to the mineral forming the walls, and diminish dissolution rates along the flow pathways.The increase of concentration within each conduit is described by a differential equation derived from a mass balance within an infinitesimal segment of conduit: where F(x) is dissolution rate at a coordinate x along a conduit, Q the flow rate and P(x), the conduit's perimeter at x.
Integration of Eq. ( 11) along a conduit gives the amount of rock dissolved within the conduit.
The dissolved load is added to the downstream junction of the conduit and is then treated as a conservative tracer by the pollutant routing code of SWMM.
In most scenarios presented in this work, transport controlled dissolution prevails.Therefore, dissolution rates are dependent on the flow velocity.A case, where the dissolution rates are almost entirely surface controlled, is also presented.

Dissolutional enlargement
Dissolution rates are rates of dissolutional enlargement v in [LT -1 ].In pressurised conduits, the cross-section changes uniformly during dissolution (Fig. 5).In a time step Δt, a conduit enlarges by vΔt, while its centre remains at the initial position.For a conduit with a free surface flow, only the wetted part of the wall is dissolved.Therefore, a transition from tube to canyon-like channel is expected.Although SWMM allows arbitrary channel geometries, the tube shape is used also during the vadose conditions in our model.To this extent an approximation is used, where the bottom of a conduit with a free surface flow incises with the true rate v and its radius increases with rate k•v, where k is the wetted fraction of the conduit perimeter.The centre of the conduit lowers with the rate (1-k)v.

The model structure
Two basic settings are presented: first a model of a Low-dip network is presented as conceptually shown in Fig. 3.This scenario is used to examine the evolution of conduit network in a plan view.In a second scenario, a highly inclined High-dip network is modelled to explore the vertical organisation of flow pathways, or evolution of the conduit network in dimension of length and depth (sensu Ford and Ewers (1978)) .
Figure 6 introduces a model structure for the Low-dip network.Circular conduits with length L and initial diameter D are assembled in an inclined rectangular grid.The orientation of the grid plane is marked geographically, N, E, S and W. All conduits are 10 m long, with initial diameters on the order of a few millimetres.Water enters the system through selected junctions indicated by arrows on Fig. 6a and flows out on the eastern boundary.Figure 6b presents junction geometry: each junction is defined by an invert elevation h 0 , relative to the base level, an inlet offset h c, which is the elevation of the conduit inlet relative to the invert elevation, and h max , the maximal depth of water in the junction.If the hydraulic head at a junction is above h max , the junction surcharges.
Figure 6c shows a side view of the model.The invert elevations increase from E to W, 1 m per junction.The slope of the W-E conduits is therefore 0.1 and N-S oriented conduits are horizontal.The inlet offset defines how much a conduit can incise.To keep conduits from bottoming out as they incise the inlet offsets, h c , are set to a large value of 100 m.Maximal depth at junctions h max is 120 m for all, except for the input junctions where h max is 111 m.
There is no storage at the junctions.
Each of the junctions on the E boundary is connected to a large conduit (D = 5 m) that freely drains water to the outfall (see Fig. 6c).These conduits play no role in the network genesis.
Their role is to effectively drain all the water arriving to the E junctions.The inverts of these junctions are at the base level and so is the inlet of the outfall conduit.This way the junctions on the E boundary allow a free outflow of the system along that face.
In the High-dip model (Figure 7), the slope of the network (and therefore the conduits) is 0.99 from top to bottom and 0.1 from left to right.We use the expressions vertical for the steep conduits and horizontal for the gradual ones.Water enters on the top side and exits at the seepage face on the right side.The bottom and left boundaries are impermeable.In all junctions, gradual (horizontal) conduits are positioned 1 m above the steep (vertical) conduits, which assures preferential flow along the vertical plane in vadose conditions (see Fig. 7b).
Flow along the horizontal conduits is active only when the junction is flooded above their inlets.The outflow is realised as in the Low-dip case, with large conduits connecting junctions to outfalls on the right boundary.

Low-dip networks
We start with a simple scenario where all conduits have the same length (10 m), the same initial diameter (0.005 m) and the same inlet offsets.The network dips from W towards the free outflow boundary on the E side with the slope 0.1 .The model is run for 50 steps of 300 s, in total 15 000 s.The rock used is salt.Nevertheless, most of the flow from upstream inputs occurs along a direct line of W-E oriented conduits, which evolve most efficiently (Fig. 8c and d).On Fig. 8e, the In3 has become vadose and in a similar manner now attracts flow from In4 and In5.However, the direct line connecting In4 to the boundary takes most of the flow and grows most efficiently.
Figure 8f shows the final stable flow configuration.All the inputs drain the available recharge, with the direct pathways between the inputs and the E boundary being the only ones that contain active flow.
A detailed look at Fig. 8 reveals that at any time, looking at the conduits draining a particular node, the highest flow rates are along W-E conduits, which consequently evolve more efficiently than other conduits.The inlet offsets of W-E conduits incise faster than others and eventually the water level at the junction falls below the lower edges of the other conduits, leaving only the W-E conduits active.This is schematically shown on Fig. 9a, where two outlets from a junction are compared; outlet 1 evolves more during the phreatic stage and, therefore, the bottom of the conduit reaches a lower elevation.Consequently, outlet 1 ultimately captures all water during the vadose entrenchment.Several other realisations of this scenario with different recharge rates at the inputs have ended with the same final distribution of active conduits.
At this point a short note is needed to explain what is meant by a stable flow configuration.In the case of constant recharge, the configuration is considered to be stable when all junctions are drained by one conduit only, i.e. there are no downstream bifurcations remaining.This is the case in Fig. 8f.In most of the other presented model runs a few outflow bifurcations remain at the last presented timestep.These bifurcations would eventually die out if the model was run long enough.We will use the term quasi-stable to describe such situations.
The next step towards less idealised scenarios is to assume that the initial inlet offsets of conduits are randomly distributed within the range of 1 m. Figure 10 shows the network when a quasi-stable flow pattern has been established, which is now more complex than in the previous case.The general evolution is similar, progressing upstream, but some N and S oriented conduits may have initial inlets low enough to keep the lowest position until the vadose transition occurs and they capture all the flow from a junction.This is schematically illustrated in Fig. 9b. Figure 11 presents the evolution of a network with initial conduit diameters drawn from a uniform distribution with a range of 10 -4 m to 10 -2 m.Initial offsets are the same for all nodes.
Generally, the evolution follows the concepts described in Fig. 8.In the pressurised phase, the selection of efficient pathways depends also on the conduit diameters and the W-E conduits are not necessarily the ones with the highest flow rates.
Figure 12 shows the evolution of total discharge from the network over time.Initially, most of the available recharge flows over the surface.First In1 and In2 integrate with full recharge summing 2 m 3 /s.After the gradient for In3 is increased, In3 integrates and the discharge rises to 3 m 3 /s.Then pathways from In4 and In5 start to contribute as these two pathways integrate.
Another selection mechanism becomes active at the transition to a free surface flow, which is shown on Fig. 13, where, a few snapshots of the SW part of the network show the evolution of several competing pathways evolving from input In5.The junctions of interest are marked by 1 to 3 and enclosed in grey circles at 4800 s.In the pressurised flow regime (4800 s), the N-S oriented conduits, marked by a, grow faster than the W-E oriented conduits marked by b at all three junctions, because conduits a belong to pathways with smaller resistance to flow.
When the flow is pressurised, the flow partitioning between two competing pathways, connecting the same junctions is divided based on the resistance to flow.Note that conduits b are parallel to the dip of the network, while conduits denoted by a are perpendicular to it.The slope of individual conduits and the distribution of slopes along the pathways plays no role.This is not the case in a free surface flow regime, where the slope of the conduit that drains the node is important.When a junction becomes vadose, the flow out of the junction through initially larger, but less steep conduits can be redistributed to more favourable steeper conduits.This leads to downstream redistribution of flow which can make part of the network inactive or change the flow from pressurised to free surface or vice versa in some of the conduits.The described situation is schematically shown on Figure 14, where two pathways, Figure 15 presents a quasi-stable flow and network pattern for the case identical to the one presented in Figure 11, but where the plane of the network is additionally tilted from N to S for 0.3 m per node.The tilting makes flow towards S preferential to flow towards N, which is clearly seen in the resulting pattern.The input In4 now joins In3.Because it is near the boundary, the input In5 has no option to develop towards S, except that the pathway heading S from the input (conduit a at In5 in Fig. 13) now persists much longer.
Other scenarios with more complex settings, such as networks with 50 x 50 nodes and networks with irregular recharge, were modelled and additionally confirmed the observations given above.
Finally we turn to a network where dissolution rate is dominantly surface controlled, as is supposed to be the case for limestone.To this end we have modelled a network, identical to the one in Fig. 11, but with α s , c eq and D set so that dissolution rates are several orders of magnitude smaller and almost entirely depend on the saturation state of the solution rather than flow velocity.Since the system is in the post-inception stage the ratio of discharge to flow length (Q/L) in many flow pathways is high enough that they evolve with the maximal growth rates.All conduits and channels along these pathways incise with the same rate.Fig. 16 shows the situation at 500 y, when a quasi-stable flow pattern has evolved and the complete network is vadose.All active channels with flow have almost the same inlet offsets and the same incision rates.Note that the colours tell the rate of increase of diameter, which is a product between dissolution rate (which is very uniform in case of surface controlled rates) and the fraction of conduit being flooded.Therefore, colours in this Figure mostly tell how full the conduits are; see also discussion in Section 2.4.The resulting flow pattern is, aside from the initial distribution of diameters and boundary conditions, a consequence of two rules: 1) at each node, channels aligned oriented with the dip drain more flow than channels perpendicular to the dip, 2) if only horizontal channels drain the node, flow is distributed evenly.The presented scenario is highly idealistic and the results and interpretation should be taken with care.In nature, the dissolution rates change with changing lithology, the initial offsets are not even, sediments can play important role, and we may question if purely surface controlled rates are reasonable.However, the model supports the ideas of Palmer (Palmer, 1991) , that maze caves develop in situations where Q/L is large along many alternative routes.

High-dip network
We now turn to the situation where the network is steep (almost vertical).As this network presents a vertical cross-section of karst, we omit the geographical notation and use top, bottom, left and right for the sides of the networks.
Similar models for laminar flow have been presented by Gabrovsek & Dreybrodt (2001) and by Kaufmann (2003).The basic result of these prior models was a continuous drop of the water table due to increased transmissivity of the network and the formation of base level conduits.If a fixed head boundary was applied, competition between a high conductivity zone along the water table and prominent conduits within the phreatic part of the network resulted in a complex pattern of evolved conduits.For many more scenarios of this modelling approach the reader is referred to the book of Dreybrodt et al. (2005)

The homogenous case with recharge distributed over the top nodes
Figure 17 presents a case where all conduits are 10 m long with initial diameter of 0.005 m.A maximum possible recharge of 5 l/s is distributed to all input nodes (blue arrows on Fig. 17a) on the top.The left column shows flow rates as line thicknesses and colours, as denoted in the legend, at five different time steps.Although the term "water table" might not be applicable for such discrete networks, we will use it for the line along the highest flooded nodes (dotted blue lines in Figs. 17 c and d).The right column shows the conduit diameters as coded in the colour bar for each figure.Equipotential lines in the left column show the distribution of hydraulic head, given in meters.
Initially (Fig. 17a), a small part of the available recharge enters the network.At the top-right all the recharge is drained directly into the outfall junction (marked by a red circle on Fig. 17a.The flow rates within the conduits are small and dominant along the vertical conduits (top to bottom).Flow along horizontal conduits is small and increases from left to right.
After 600 s (Fig. 17 b) the entire network is still pressurised.Horizontal conduits have evolved sufficiently to drain more flow brought in by initially developed vertical conduits.
Accordingly, the potential gradient becomes oriented to the right and is highest close to the boundary.Conduits at the top-right corner experience fastest growth and capture almost all recharge from the inputs.The flow in the left part of the network is small and the hydraulic potential field is relatively flat there.After 1200 s (Fig. 17 c) the top-right corner has become vadose).In this area, the recharge is carried vertically to the water table.The flow rates are highest along the water table and diminish with distance from it.
However, widening is still substantial below the water table which additionally increases the network permeability and downwards retreat of WT.The process continues until the WT drops to the base level and only vertical recharge conduits and a master conduit at the base continue to grow.The vertical conduits have been widened through the entire evolution, the uppermost for the longest time and they are therefore largest.The diameters decrease from top to bottom.On the other hand, the diameter of horizontal channels increase from left to right, as they evolve only below the water table.Therefore, deeper conduits have more time to evolve.

Inhomogeneous case
In the case shown on The initial diameter of the top horizontal line of conduits is 0.1 m.
A recharge of 100 l/s is introduced to the top-left junction (see the blue arrow on Fig. 18a. The two given legends for flow rates and diameters are valid for all figures.At 3000 s (Fig. 18a), about one fourth of the available recharge is captured and drained directly to the outfall by the top line of horizontal conduits.
Pathways along the conduits with initially larger diameters evolve efficiently and capture an increasing amount of flow.
At 9000 seconds (Fig. 18b) about 70% of the flow is captured by the junction marked by a blue triangle and denoted by 1 in Fig. 18b.It feeds a line of vertical conduit that discharge into outflows through horizontal conduits.Numbers on the conduits in the top-right region denote flow along the conduits in l/s.The discharge to the outflow diminishes downwards.
However, these conduits widen effectively and cannot sustain a pressurised regime, so that the position of highest outflow migrates downwards.
By 24000 s, the outflow position has retreated to the bottom (Fig. 18c).When the vertical pathway downwards from point 1 becomes vadose, it provides a free outflow boundary and triggers the development of pathways draining sink points 2 and 3 (Fig. 18b,c), which soon capture all the flow.On Fig. 18c, the flow along the top line has retreated to point 3 and throughout the remainder of the simulation continues to retreat towards the left to points 4 and 5 (Fig. 18d) .Ultimately, the flow is captured by the node at point 5 (Fig. 18e).Similarly, the flow migrates from top to bottom, towards the deeper connecting pathways.Figure 18e shows the stable flow situation at 75000 s, where all the flow follows one single pathway.
Downward and leftward progress is slow because some of the conduits to the left are initially small and the permeability is low.In comparison with a uniform network with distributed recharge, the development follows initially prominent pathways, with progressive upstream flow capturing.Soon after a pathway becomes vadose, the flow is overtaken by the evolving pathways to its left.

The role of prominent structures
The progression mechanism described above, is demonstrated clearly by a final idealised, but telling, example.We assume three vertical conduits ("wells") with an initial diameter of 0.2 m, extending completely through the domain in the vertical direction.
These are connected with 5 evenly spaced horizontal conduits with initial diameter 0.005 m extending across the domain.All other conduits are effectively impermeable, with a diameter of 10 -5 m.A maximum possible recharge of 100 l/s is available to the prominent vertical conduits (wells) as marked by the arrows at the top of Figure 19a.
Initially (Fig. 19a), all conduits are pressurised.There is almost no gradient left of W3, where evolution is slow or none.High gradients exist between W3 and the outfalls, the highest being along the deepest horizontal conduit, which has the highest flow and evolves most efficiently.
As W3 becomes vadose, it presents a free outflow boundary for the flow from its left and the gradient along the horizontal conduits connecting W2 to W3 builds up.These conduits now experience fast evolution with rates increasing from the top to the bottom (Fig. 19b).The mechanism progresses leftwards: when W2 becomes vadose, W1 connects to it as shown in In the pressurised phase, the flow out from a junction is distributed to the outlet conduits, according to their resistance to flow and the distribution of hydraulic heads.This also defines the rate of their inlet incision.When a junction becomes vadose, the conduit with the lowest inlet entry elevation has an advantage and is a candidate to take all the flow.However, under vadose conditions the conduit's alignment with respect to the dip of the network becomes important, as higher slope generally invokes higher energy grade, higher flow velocity and faster incision.A conduit that gains advantage in pressurised conditions, can be surpassed by a conduit with a higher slope, which has an advantage in free surface conditions.Once the stable flow pattern is established, the flow follows a system of conduits that all occupy the lowest position in their upstream junctions.

High-dip scenario
In a homogenous scenario, the evolution is focused to the transitional area between pressurised and free surface flow, the "water table".The flow from the surface is gravitational along the vadose channels down to the water table.There, it is largely focused to the conduits close to the water table.The scenario demonstrates a relatively smooth drawdown of the water table due to increasing permeability in the phreatic zone.The end result is a relatively uniform network with a growing base level conduit.Similar results were obtained by Gabrovšek and Dreybrodt(2001) and by Kaufmann (2003), where only dissolution in the phreatic zone was considered.
The inhomogeneous case demonstrates the evolution of invasion vadose caves based on flow diversion.The drawdown of the phreatic zone is irregular, following fast evolution of prominent pathways and progressive upstream flow capturing.Such a scenario can produce extended an network of steep vadose passages.
Deeply penetrating conductive structures can play an important role as they transfer surface water deep into the massif and redistribute hydraulic gradients.This way fast evolution along deep horizons can be triggered.

Conclusion
The On the other hand the model opens new challenges related to evolution of karst aquifers in vadose settings.Further work is needed to improve estimation of dissolution rates and the related role of sediment transport and mechanical erosion.Further steps towards more realistic modelling domain and boundary conditions are also needed.In fact, a single low-dip plane is a scenario, which is not common in the nature.A careful step towards 3D models that simulate speleogenesis in both, phreatic and vadose conditions is therefore needed.By careful, we mean gradual adding of complexity, so that at each new step all mechanisms from previous steps are well understood.The presented model allows such extensions.
At the same time, we have to keep in mind the modelling results are not stand-alone, they should progress hand in hand with conceptual models based on the field observations.Grey lines show scenario with constant input at 1.5T B (Fig. 1e).

2.
Calculate flow in the network.3. Couple flow, dissolution and transport to calculate dissolution rates in all conduits.4. Change the conduit diameter within a time step according to the dissolution rate and return back to Step 2 or exit the loop when a stable flow pattern is established or no substantial changes in flow pattern are expected.

Figure 8
Figure 8 presents six snapshots of the network's evolution.Five inputs with Q max = 1000 l /s are marked by circles and denoted by 1-5 on Fig. 8a.The left column shows flow rates and flow directions.Flow rates are denoted by line thicknesses and flow directions by colour; red represents flow towards N or W and black towards S or E. If the flow is pressurised, the colours are saturated; pale colours denote conduits with free surface flow.The right column and b connect two nodes.Pathway a is initially larger, drains more flow, and widens more efficiently in the pressurised phase.When the conduit turns vadose, the flow rates in a drop due to the low slope of the channel as it leaves the junction.If, at the transition to free surface flow, the water level in the upstream node has not dropped below the inlet of pathway b, the steeper entry into pathway b as it leaves the junction causes b to incise faster and progressively capture more flow.
Figure 18  we assign a more complex distribution of initial conduit diameters.The initial diameter (d o ) of each conduit is constructed as a sum of a group contribution (d g ) which is given to all conduits aligned along the same line, and an individual contribution (d i ).These are both random, sampled from a uniform distribution, where that conduits along a certain line get the individual contributions is 0.5.Using this group contribution, we enhance the potential importance of conductive structural lines.

Fig 19c .
Fig 19c.In Figure 19d, a stable flow condition is shown, where all the flow follows the wells which feed the base level channel.
presented model closes some of the open questions, which have not been addressed by the older existing models.The final flow pattern results from all stages of network development, starting with the initial stage, continuing with the growth, integration and expansion under pressurised flow and, what is demonstrated by this model, with the final selection of stable flow pathways on a local scale during and after transition to free surface flow regime.

Figure 1 :
Figure 1: Evolution of 2D fracture network under pressurised flow.Panels show aperture widths and dissolution rates at different stage of evolution.Size of the domain is 1 km x 1 km, initial aperture width a 0 = 0.02 cm, except for the line P1 , where a 0 = 0.03 cm.Linear and forth order dissolution kinetics for the limestone is used (see Dreybrodt et al. (2005) for details).

Figure 2 :
Figure 2: Profile of hydraulic head along pathways P1 (dashed lines) and P2 (full lines) from Fig.1.Profiles are taken at different time steps, given in units of breakthrough time (T B ).

Figure 3 :
Figure 3: Conceptual framework.A conduit network with point recharge at selected locations indicated by arrows.Recharge is limited by the position of the land surface h max or by maximal available recharge Q max .

Figure 4 :
Figure 4: The use of a Preissmann slot enables use of the same set of equations for conduits with free surface flow and conduits with pressurised flow.

Figure 5 :
Figure 5: Growth of a conduit with pressurised flow and a conduit with free surface flow.r is radius, k is the fraction of wetted perimeter, v incision (growth) rate.

Figure 6 :
Figure 6: The model structure for the Low-dip network.a) A conduit network with discrete water inputs, marked by arrows.Boundaries are denoted geographically.Outputs are along the E boundary.b) Geometry and parameters of a junction.c) The side view of the model, also showing a large conduit connecting E junctions to an outfall.

Figure 7 :
Figure 7: The model structure for the High-dip scenario.a) The slope of the network is 0.99 in from top to bottom and 0.1 from left to right.The right boundary is a seepage face with free outflow.Inputs are on the top.b) Junction geometry: high-dip ("vertical") conduits are positioned below the low-dip ("horizontal") conduits.

Figure 8 :
Figure 8: Six snapshots of the evolution of Low-dip network with uniform initial diameters and inlet offsets.Left: flow rates (width) and flow direction (Red = flow towards E or towards N, Black/Grey = flow towards W or towards S).Right: diameters (width) and widening rates (colour).The codes below show thicknesses, flow rates and widening rate.The values at the bar codes correspond to the thickest lines in the flow rate and diameter bars and to the warmest colour in the bar for the widening rate.The scales are linear with the thinnest lines and dark blue colours representing no flow, no widening the and smallest initial diameter.

Figure 9 :
Figure 9: Left: the geometry of a junction.Right (a and b): Scheme of two outflows during pressurised flow (top) and free surface flow (bottom).a) initial inlet offsets for both outflows are equal.b) Initial inlet offset of outflow 2 is smaller so that the outflow has a lower elevation.Blue arrows indicate the amount of flow drained by each outflow, and the blue shading indicates the water table.

Figure 10 :
Figure 10: A network with uniform initial diameters and initial inlet offsets randomly distributed within vertical span of 1 m.

Figure 11 :
Figure 11: Evolution of a Low-dip network with randomly distributed initial diameters.

Figure 12 :Figure 13 :
Figure 12: The time evolution of total discharge from the network in Fig. 11.

Figure 14 :
Figure 14: Distribution of flow between two pathways depends on the flow resistance when the flow is pressurised.The pathway a has with lower flow resistance grows faster.After the transition to free surface flow, the pathway b with higher exit slope from the junction can capture more flow and incise faster.

Figure 15 :
Figure 15: Quasi-stable state of network with same structure as presented in Fig. 11, but the plane of the network is additionally tilted from N to south, for 0.3 m per node

Figure 17 :
Figure 17: Evolution of homogenous sub-vertical network.Blue arrows on Fig. 17a denote inputs.Isolines and values present the hydraulic potential [m]

Figure 18 :
Figure 18: High-dip network with random initial distribution of conduit diameters.Flow enters at the top-left edge of the network as pointed by a blue arrow.Values on Fig. 18b show flow rates along the selected individual conduits

Figure 19 :
Figure 19: High-dip network with three prominent conduits (wells), marked by W1 to W3.A recharge of 100 l/s is available to the prominent conduits.

Table 1 :
List of rate constants and other parameters used in this work.
* To have dissolution rates expressed as a velocity of wall retreat, concentration [NL -1 ] is multiplied with molar mass [MN -1 ] and divided by the density [ML -3 ] of the mineral forming the rock and being dissolved.This makes c eq it dimensionless.List of Figures: