FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES ODVISNOST KRAŠKEGA RAZVOJA RAZPOKE OD VRSTE IN ZAPOREDJA VODOTOPNIH KAMNIN, SKOZI KATERE POTEKA Georg KAUFMANN1*, Franci GABROVŠEK2 & Douchko ROMANOV3* Izvleček UDK 551.435.8 Georg Kaufmann, Franci Gabrovšek & Douchko Romanov: Odvisnost kraškega razvoja razpoke od vrste in zaporedja vodotopnih kamnin, skozi katere poteka Z numeričnim modelom smo raziskovali razvoj sekundarne poroznosti v enodimenzionalni razpoki v apnencu, sadri in an- hidridu. Cilji raziskav so bili: 1) ugotoviti razlike med razvojem plitvih in globokih razpok v različnih kamninah; 2) oceniti učinek morebitnega izločanja enega od mineralov v razpoki in 3) dinamika rasti razpoke, ki gre skozi plasti različnih kamnin. Rezultati kažejo na primerljiv razvoj razpok v apnencu in sadri, pri čemer je čas razvoja v obeh kamninah precej različen. V anhidridu je zaradi drugačne kinetike raztapljanja razvoj še hitrejši. Izločanje ob spreminjajočih se geokemičnih pogojih lahko hitro zamaši razpoke in posebej v večjih globinah spre- meni poti razvoja mreže kraških prevodnikov. V primeru, ko so razpoke v apnencu prekrite s plastjo sadre, kar v naravi pogosto opažamo, je hitrost širjenja najhitrejša, kar daje vodnim potem vzdolž takih razpok izrazito primerjalno prednost pri zajema- nju razpoložljivega toka. Ključne besede: vodotopne kamnine, enostavna razpoka, raz- tapljanje in izločanje, incepcijski horizont. 1 Institute of Geological Sciences, Geophysics Section, Freie Universität Berlin, Malteserstr. 74-100, Haus D, 12249 Berlin, Germany, e-mail: georg.kaufmann@fu-berlin.de 2 Karst Research Institute ZRC SAZU, Postojna, Slovenia, e-mail: gabrovsek@zrc-sazu.si 3 Institute of Geological Sciences, Geophysics Section, Freie Universität Berlin, Malteserstr. 74-100, Haus D, 12249 Berlin, Germany, e-mail: douchko.romanov@fu-berlin.de *Corresponding Author Received/Prejeto: 10.05.2017 COBISS: 1.01 ACTA CARSOLOGICA 46/2−3, 199–216, POSTOJNA 2017 Abstract UDC 551.435.8 Georg Kaufmann, Franci Gabrovšek & Douchko Romanov: Fracture evolution in soluble rocks: From single-material fractures towards multi-material fractures We employ a numerical model describing the evolution of sec- ondary porosity in a single fracture embedded in soluble rock (limestone, gypsum, and anhydrite) to study the evolution of isolated fractures in different rock types. Our main focus is three-fold: The identification of shallow versus deep flow paths and their evolution for different rock types; the effect of precip- itation of the dissolved material in the fracture; and finally the complication of fracture enlargement in fractures composed of several different soluble materials. Our results show that the evolution of fractures composed of limestone and gypsum is comparable, but the evolution time scale is drastically differ- ent. For anhydrite, owing to its difference from calcite in the kinetical rate law describing the removal of soluble rock, the evolution is even faster. Precipitation of the dissolved rock due to changes in the hydrochemical conditions can clog fractures fairly fast, thus changing the pattern of preferential pathways in the soluble aquifer, especially with depth. Finally, limestone fractures coated with gypsum, as occasionally observed in caves, will result in a substantial increase in fracture enlarge- ment with time, thus giving these fractures a hydraulic advan- tage over pure limestone fractures in their competition for cap- turing flow. Key words: soluble rock, single fracture, dissolution and pre- cipitation, inception horizon. ACTA CARSOLOGICA 46/2–3 – 2017200 GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV Soluble rocks such as limestone, dolostone, gypsum, and anhydrite are characterised by the dissolution of material by water, often enriched with carbon dioxide, sulphides or other organic acids. Removal occurs both on the rock surface and along fractures and faults in the subsurface in the case of a telogenetic evolution of the rock. While in insoluble rocks the permeability given by the continuous subsurface flow paths between interconnected fractures and bedding planes remains mainly constant, in soluble fractured aquifers permeability increases substantially with time. The resulting flow paths in soluble fractured aquifers provide efficient drainage paths for water, and surface flow often completely disappears, as the enlarged fractures can accommodate large amounts of water, even during larger recharge events. The subsurface flow paths are, however, preferentially developed, and often voids enlarged to the meter size and more provide access to caves carrying entire underground rivers. We focus on chemical reactions during the evolu- tion of a fracture in soluble rocks by comparing different rock types (limestone, gypsum, anhydrite). For this pur- pose, we simplify our fracture to a single, isolated circu- lar conduit with variable diameter, through which flow is driven either under laminar or turbulent flow condi- tions, with a conduit roughness coefficient mimicking small-scale wall irregularities, and transport of dissolved species as an advective process. The interaction of this single isolated conduit with other conduits in the aquifer is therefore neglected. Szymczak et al. (2006; 2009; 2011) have emphasized that a single fracture modelled in two dimensions will lead to a break-up of the dissolution front, inducing fin- ger flow and a preferential enlargement along the two- dimensional fracture. Thus, the two-dimensional frac- ture will evolve more rapidly than its one-dimensional analogon. However, this behaviour, termed exchange flow between fractures in a karst aquifer, has been stud- ied extensively in the past, and we will argue that the ef- fect of exchange flow will be similar in all chosen rock types and thus is already understood. In this paper, we address three key points: (i) Dis- cuss shallow and deep flow and evolution in fractures for different soluble rock types (limestone, gypsum, an- hydrite), (ii) include precipitation into the evolution of the fractures, (iii) discuss the evolution of fractures com- posed of several soluble rock layers. We have organised the paper as follows: In section 2, we discuss processes controlling the evolution of con- duits, which we want to address by our numerical mod- elling. In section 3, we introduce physical and chemi- cal principles for flow and evolution of a single conduit embedded in soluble rock. In section 4, we discuss the modelling results, starting with the enlargement of sin- gle conduits in limestone, gypsum and anhydrite in shal- low and deep conditions. We then proceed to describe changes due to possible precipitation of the soluble rock. We finally look at conduits composed of several soluble rock types. In section 5, we discuss and summarise our results. INTRODUCTION PROCESSES In this section, we discuss mechanisms for fracture widen- ing, fracture clogging, and the interaction of different dis- solved species in fractures. Fracture enlargement in soluble rocks is not a uniform process, but highly selective. The evolution of caves with respect to the water table was a debated phenomena in the mid-20th century, largely made upon hydrogeological principles. Ford and Ewers (1978) have reconciled these principles with the definition of the four-state model: In this model, fracture spacing controls the type of cave evolution. For low frac- ture spacing, caves with deep phreatic loops can evolve, higher fracture spacing facilitates evolution of caves along the water table. The four-state model is capable of explaining cave geometries for different geological settings, such as the evolution of caves along the water table, if fracture spacing is dense, and the evolution of deep bathy-phreatic cave loops, if the fracture spacing is less dense. Worthington (2001) came up with a different expla- nation of deep phreatic flow paths: Viscosity decreases with increasing temperature, thus with depth, and there- fore flow deeper in the aquifer can have a slight hydraulic advantage over shallow flow (hydraulic control). Kauf- mann et al. (2014), however, have shown by means of numerical modelling that the flow enhancement due to the decrease in viscosity with depth (hydraulic control) is counteracted by the reduction of solubility with depth (chemical control) and the reduction of initial fracture width with depth due to lithostatic pressure (structural control). ACTA CARSOLOGICA 46/2–3 – 2017 201 FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES MODELLING FLOW AND EVOLUTION The evolution of fractures in soluble rock is controlled by both flow through the fracture and changes of the geom- etry of the fracture due to dissolution of and precipita- tion on the fracture walls. These two groups of proper- ties, which we term hydraulic properties and chemical properties, are coupled, and this coupling often results in a feedback mechanism responsible for the preferential enlargement of fractures and bedding partings in soluble rock. In this section, we collect the relevant relations needed to describe and understand this feedback be- tween hydraulic and chemical processes. We will use the mathematical concept of a circular conduit to mimic the real fracture. In Fig. 1, a typical cross section 3 km long and 2 km deep for such a conduit is outlined: Two conduit paths are shown (white with black outline): one flow path close to the surface along the water table, the other flow path resembling bathyphreatic conditions and reaching into zones, where properties such as tempera- ture, pressure and related parameters differ substantially from their surface values. HyDRAULIC PROPERTIES We start by defining a circular conduit of diameter d [m] and length l [m], and the deepest point of the conduit as zmax [m]. Flow through this conduit can be described in general by the Navier-Stokes equations, which in our case can be simplified to classical Darcy flow under laminar or turbulent conditions (e.g. Beek & Muttzall 1975): Recently, the hypothesis of inception horizons (e.g. Lowe 1992, 2000; Filipponi & Jeannin 2006; Filipponi et al. 2009) extended the idea of preferential enlargement of fractures and bedding planes. In this hypothesis, an inception horizon is identified as part of the soluble rock, which has different physical, lithological or chemical properties than the surrounding rock. Often, inception horizons identified in the field are covered with gypsum. We will explore the evolution of a single conduit in lime- stone coated with gypsum to explore the importance of chemical control in conduit evolution, beside the hy- draulic control described by the four-state model. The clogging of fractures can also play a key role in the development of flow through aquifers. Beside clogging through sediments, large lithostatic pressure, or bacterial activity, soluble rocks can in places deliver solution supersaturated with a mineral species, which then starts to precipitate and reduce fracture width and thus flow. As an example, we report the failed geother- mal drill holes in Staufen (Breisgau, Germany), where the conversion of anhydrite to gypsum, once the anhy- drite came into contact with water, caused widespread damage. In 2007, seven boreholes for geothermal heat exchange were drilled in the city of Staufen to a depth of 140 m (e.g. Sass & Burbaum 2010; Lubitz et al. 2013). One borehole connected two independent groundwater horizons and flow through an intermittent anhydrite lens caused widespread dissolution of the anhydrite and subsequent precipitation of gypsum. The larger volume of the precipitated gypsum caused surface uplift of up to 1 cm/month in an area 100×150 m, which caused wide- spread structural damage to the historical buildings. We will explore the interaction between anhydrite-gypsum in a single fracture. fig. 1: Model setup for the karst conduit evolu- tion. temperature is shown as color contours, ranging from 10 °C at the surface to around 60 °C in 2 km depth. The corresponding prop- erties of water are shown as black lines: density ρ (dashed line temperature effect only; solid line temperature and pressure effect), viscosity η, pressure p (dashed line water pressure, solid line lithostatic pressure), calcium equilibrium concentration ceq (dashed line temperature effect only, solid line temperature and pres- sure effect; not distinguishable), fracture-depth relation d/di. The thick white lines with black outline indicate two conduit paths, one along the water table, one deep bathyphreatic cave. ACTA CARSOLOGICA 46/2–3 – 2017202 with Q [m3/s] the flow rate through the fracture, ∆h [m] the drop in hydraulic head between both ends of the fracture, g [m/s2] gravitational acceleration, and f [–] the friction factor. Note that (1) is a reformulation of the classical Poiseuille law for both laminar and turbulent conditions. The friction factor can be defined as (e.g. Jeppson 1976): with Re [–] the Reynolds number, and w [m] the wall roughness, a parameter describing the roughness of the conduit walls. The different representations of f de- scribe laminar flow (fl) and different stages of turbulent flow (fs-smooth turbulence, ft-transitional turbulence, fr-rough turbulence). For the latter, the largest value con- trols turbulent flow behaviour. Note that both Reynolds number and friction factor have to be calculated with a standard iterative procedure. The friction factor depends on the dimensionless Reynolds number Re: with A [m2] the cross-sectional area of the conduit, ρw [kg/m 3] the fluid density, and η [Pa s] fluid viscos- ity. If the Reynolds number is below a certain threshold (Rec ~ 2200), flow in the conduit is laminar, above the threshold flow becomes turbulent. Thus, depending on the flow regime, we arrive at a linear relation between hydraulic head drop and flow, describing laminar flow, and a quadratic relation, repre- senting turbulent flow conditions. The flowrate depends strongly on the conduit width (power-law of the order 3−4). The change in conduit width can be described (e.g. Kaufmann et al. 2014): with ti and ti+1 [s] two consecutive time steps, f [mol/m 2/s] the calcium flux rate, mrock [kg/mol] the atomic mass of the dissolved rock, ρrock [kg/m 3] the density of the dis- solved rock, and ∆t = ti+1−ti [s]. Note that d(t = 0) = d0 is the initial width of the conduit. The change in conduit width therefore depends on the calcium flux rate, with f > 0 describing enlargement of the conduit size by re- moving material from the conduit walls, and f < 0 reduc- tion of the conduit width by precipitating material on the conduit walls. The calcium flux rate itself is a function of the calcium concentration c [mol/m3] in the conduit. During dissolution, c will increase along the conduit be- cause of the removal of soluble rock from the conduit walls, until c reaches the calcium equilibrium concentra- tion ceq [mol/m 3]. The change in calcium concentration along the con- duit is calculated as (e.g. Dreybrodt et al., 2005) with xi and xi+1[m] two neighbouring points along the conduit, ∆x = xi+1−xi, f(x) [mol/m 2/s] the calcium flux rate, and P(x) [m] the perimeter of the conduit. Note that c(x0)=cin is the input calcium concentration cin [mol/m 3] entering the conduit. CHEMICAL PROPERTIES The dissolution or precipitation of material from a con- duit in soluble rock depends on the amount of soluble material that can be dissolved from the rock wall into the solution flowing through the conduit. Depending on the type of soluble rock, the equilibrium concentration of the dissolved species is a function of several properties. We will focus on gypsum (CaSO4 · 2H2O), anhy- drite (CaSO4), and limestone (CaCO3) as soluble rock types. The chemical reactions for these different soluble rocks in an aqueous solution can be described through the following set of reactions (e.g. Duan & Li 2008; Li & Duan 2011): The first reaction describes the dissociation of water, the second the solution of carbon dioxide (CO2) in water, the next two ones the dissociation of carbonic acid, and the remaining equations the solution of gypsum, anhydrite and limestone, respectively. All of the above equations can be described as equi- librium reactions with their respective mass action co- efficients, ki. The calcium equilibrium concentration for the soluble rock types can be derived as approximate function to excellent accuracy (e.g. Dreybrodt 1988): GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV ACTA CARSOLOGICA 46/2–3 – 2017 203 with kA(t,z) (Blount & Dickson 1973) the equilibrium constant for the dissolution of anhydrite, kg(t,z) (Blount & Dickson 1973) the equilibrium constant for the disso- lution of gypsum, kh(t,z) (Weiss 1974) the equilibrium constant for the dissolution of atmospheric carbon diox- ide into water (Henry’s law constant), k0(t,z) (Wissbrun et al. 1954) the equilibrium constant for the reaction of water and carbon dioxide to carbonic acid, k1(t,z) and k2(t,z) (Millero 1979; Mehrbach et al. 1973) the equilib- rium constants for the dissociation of carbonic acid into bicarbonate, carbonate, and hydrogen, kC(t,z) (Mucci 1983) the equilibrium constant for dissolved calcite, γCa2+, γSO42− and γhCO3− the activity coefficients for calcium, sulfate and bicarbonate, pCO2 [atm] the carbon-dioxide partial pressure, and t [°C] the temperature of the solution. Our analytical expressions compare well with the calculation of the calcium equilibrium concentration derived for the full electro-neutrality condition, which we verified by comparing our ceq to results obtained from PHREEqC (Parkhurst & Appelo 2013). For limestone, (7) is valid for the open system, in which the solution is in contact with the atmosphere, and carbon dioxide will be replenished by further solu- tion of CO2 from the atmosphere. However, most con- duit enlargement takes place under closed-system condi- tions, and here the CO2 is consumed and thus decreases with dissolution. In this case, the carbon- dioxide partial pressure is (Dreybrodt 1988): with patm [atm] the initial carbon-dioxide partial pressure obtained in the atmosphere and the soil. The calcium equilibrium concentrations for lime- stone, gypsum, and anhydrite are shown in Fig. 2 as func- tions of temperature t [°C], and carbon- dioxide pres- sure pCO2 [atm] for the case of limestone. Note the large difference in ceq between limestone ceq~1−5 mol/m 3, gyp- sum ceq~15 mol/m 3, and anhydrite ceq~40 mol/m 3. Also note that ceq for both limestone and anhydrite is a ret- rograde function of temperature, while for gypsum the temperature relation is prograde for temperatures below 30 °C, retrograde above that temperature (see Fig. 2). The calcium flux rate describes flux of dissolved spe- cies from and to rock surface per unit area and per time. It is controlled by several potentially rate-limiting pro- cesses on the bedrock surface, e.g. the surface reaction at the mineral surface and the transport of the dissolved species in the solution. Flux rates have been measured experimentally for limestone (e.g. Plummer et al. 1978; Svensson & Dreybrodt 1992; Eisenlohr et al. 1999), and for gypsum and anhydrite (e.g. James & Lupton 1978; Gobran & Miyamoto 1985; Lebedev & Lekhov 1990; Je- schke et al. 2001; Jeschke 2002), and for limestone have fig. 2: Calcium equilibrium concentration as a function of temperature. Shown are curves for limestone (CaCO3), gyp- sum (CaSO4·h2O), anhydrite (CaSO4), and halite (NaCl). for limestone, the calcium equilib- rium concentration is shown for three values of partial carbon-di- oxide pressures and for open- and closed-system conditions. FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES ACTA CARSOLOGICA 46/2–3 – 2017204 been predicted numerically (e.g. Buhmann & Drey- brodt 1985a,b; Dreybrodt & Kaufmann 2007; Kaufmann et al. 2010). Flux rates f [mol/m2/s] can be described as a piece-wise function of the calcium concentration c with respect to the calcium equilibrium concentration ceq (e.g. Palmer 1991) with ki [mol/m 2/s] a rate coefficient, c [mol/m3] the actual calcium concentration, ceq [mol/m 3] the calcium equilib- rium concentration, and ni [-] a power-law exponent. The rate coefficient depends on both surface-controlled rates and diffusion-controlled rates: with ki S [mol/m2/s] the surface-controlled rate coefficient, and with ki D = 2Dceq/d [mol/m 2/s] the diffusion-control- led rate coefficient, and D [m2/s] the diffusion coefficient. For a summary of coefficients see Tab. 1. tab. 1: Parameter values for soluble rock chemistry. Limestone Atomic mass mRock [kg/mol] 0.100 Density ρRock [kg/m 3] 2600 Rate constant1 k1 [mol/m 2/s] 4x10–7 Rate constant k2=k1(1-cs) n1-n2 Rate exponent1 n1 [-] 1 Rate exponent1 n2 [-] 4 Switch1 cs [-] 0.90 Anhydrite Atomic mass mRock [kg/mol] 0.136 Density ρRock [kg/m 3] 2900 Rate constant2 k1 [mol/m2/s] 4x10 –5 Rate exponent2 n1 [-] 5.4 Gypsum Atomic mass mRock [kg/mol] 0.172 Density ρRock [kg/m 3] 2200 Rate constant3 k1 [mol/m 2/s] 1.1x10–3 Rate constant k2=k1(1-cs) n1-n2 Rate exponent3 n1 [-] 1 Rate exponent3 n2 [-] 4.5 Switch1 cs [-] 0.95 1 Buhmann et al. (1985) 2 Jeschke (2002) 3 Jeschke et al. (2001) The calcium flux rates for limestone, gypsum, and anhydrite are shown in Fig. 3 as a function of calcium concentration c for different temperatures T and carbon- dioxide pressures pCO2 in the case of limestone. For both limestone and gypsum, the flux rates for dissolution are characterised by a linear decrease with increasing cal- cium concentration, until a material-dependent thresh- old (cs) is reached. From then on, the flux rates decrease following a power law, thus are much slower (Svensson & Dreybrodt 1992; Eisenlohr et al. 1999; Buhmann & fig. 3: Calcium flux rates as a function of calcium concentration for different soluble rocks. Curves are shown for two different temperatures and in the case of limestone also for three different partial carbon-dioxide pressures. top: limestone, middle: gyp- sum, bottom: anhydrite. GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV ACTA CARSOLOGICA 46/2–3 – 2017 205 Dreybrodt 1985a,b). The reason for these non-linear flux rates at high calcium concentrations is the accumulation of impurities, which originate from insoluble material in the soluble host rock (e.g. clay). Note that for anhydrite, the flux rate is non-linear over the entire range for dis- solution, which results in drastically different behaviour, as we will see later. Once the calcium concentration c passes the cal- cium equilibrium concentration ceq, precipitation starts. Here, experimental data are scarce, and only the precipi- tation rates for limestone are based on laboratory work (Buhmann & Dreybrodt 1985b; Dreybrodt & Buhmann 1991; Dreybrodt et al. 1997). For precipitation rates of anhydrite and gypsum, we assumed a linear relation, but we also discuss potential non-linearities later. The lin- ear constant for precipitation is the negative of k1 for the specified rock. DEPTH DEPENDENCE The material parameters density, gravity, and viscosity are functions of several variables, e.g. temperature and water pressure, which will be parameterised as depth z [m]. While we can safely neglect the depth dependence of the gravitational acceleration over the depth range of even deep karst aquifers, temperature, hydrostatic and lithostatic pressure, density, and viscosity have to be pa- rameterised according to the increase in depth: with t0 [°C] the annually-averaged surface temperature, (dT/dz)geotherm [°C/m] the geothermal gradient, with ρw [kg/m 3] the density of water, and ρl [kg/m 3] the density of rock, ρ0 [kg/m 3] the respective surface values, g [m/s2] the gravitational attraction, α [1/K] the thermal expansiv- ity of water (e.g. Jones & Schoonover, 2002), κt [1/Pa] the compressibility of water (e.g. Jones & Schoonover 2002), η [Pa s] the viscosity of water, η0 [Pa s] the reference value for the viscosity of water at surface pressure and 25 °C, ai polynomial coefficients for the viscosity (Kestin et al. 1978). Note that of course the water pressure depends on the induced flow, but for the depth-dependence of pa- rameter values we use the hydrostatic pressure as simpli- fication, which is justified, as induced pressure is smaller than increase in hydrostatic pressure. A more detailed discussion of the depth-dependent properties can be found in Kaufmann et al. (2014). MODEL IMPLEMENTATION The equations developed above have been implemented by the authors into a simple numerical code in Fortran90. The program performs the following steps: Parameter values are read (e.g. d0, l, zmax, ∆h, cin, pCO2, t0, (dT /dz)geotherm). The conduit is discretised into nx-elements, and the con-• duit path is assigned, following z = zmax[1−((2x/l)−1) 2], assigning a parabolic flow path for depth. Note that spatial discretisation for the conduit is controlled by the change of calcium concentration along the conduit (see Dreybrodt et al. (2005), for a discussion of con- centration increments). A time stepping routine then calculates the evolution • of the conduit. For each time step, a parcel of water enters the first conduit element with the pre-defined calcium input concentration cin. Then for each sub-se- quent conduit element, temperature, density, viscosity, and conduit width are calculated for the given depth according to (11). Time steps are chosen small enough to ensure convergence of results. The flow rate in the conduit element is calculated, • based on (1) and depending on the flow regime, using the equivalent resistance formula for conduit elements tab. 2: Parameter values for depth dependences of material prop- erties. Geothermal gradient (dT/dz)geotherm 25°C/km Thermal expansivity α Compressibility κτ Viscosity η0 1.002x10 –3 Pa s a0 1.2378 a1 1.303 x 10 –3 a2 3.060 x 10 –6 a3 2.550 x 10 –8 tab. 3: Parameter values for standard model. Conduit length l 2000 m Conduit radius d 0.5 mm Max. conduit depth zmax 0 m Head drop Δh 10 m Surface temperature T0 10°C CO2-Pressure pCO2 0.05 atm Calcium input concentration cin 0 mol/m 3 Wall roughness w 0.00002 m Diffusion coefficient D 10–9 m2/s FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES ACTA CARSOLOGICA 46/2–3 – 2017206 in series (see Groves & Howard (1994a) and Howard and Groves (1995), for a discussion of the equivalent resistance). With these parameters, the calcium concentration (5), • the calcium equilibrium concentration ceq (7) and the calcium flux rate F (9) are calculated in each conduit element. Based on the flux rate and the time step • ∆t, the new width of the conduit element is calculated according to (4). Finally, the calcium concentration is changed accord-• ing to (5), and then passed to the next conduit ele- ment. We now have assembled all relevant information on flow, evolution, and depth dependence of the hydraulic and chemical parameter values. RESULTS & DISCUSSION In this section, we discuss results from numerical so- lutions of the flow and evolution equations for a single circular fracture. We start with a shallow conduit resem- bling a water-table cave and three different soluble rock types, limestone, gypsum, and anhydrite. We then allow the conduit to penetrate deeper into the aquifer, resem- bling a bathy-phreatic flow path. Then we proceed with a combined dissolution-precipitation of material in these different soluble rocks. Finally, we describe a limestone conduit coated with gypsum to discuss the inception hy- pothesis. DISSOLUTION IN SHALLOW SETTINGS We first discuss the evolution of a single conduit in dif- ferent soluble rocks and for three soluble rock types: This first case we term water-table conduit, as it is horizontal and close to the surface. The conduit considered has a length of l=2000 m and an initial width of d0=0.5 mm for the water-table conduit. Flow is driven from left to right, with a hydraulic head drop of ∆h=10 m, and the calcium concentration of solution entering the conduit is cin=0 mol/m 3. As climatic variables, we assume a tem- perature of T = 10 °C and a carbon-dioxide pressure of pCO2=0.05 atm. Limestone The temporal evolution of a conduit in limestone is shown in Fig. 4 (top). For the given temperature and carbon-di- oxide pressure, the calcium equilibrium concentration is around ceq 2.11 mol/m 3. The conduit evolves as a typical conduit in soluble rocks often discussed in the literature as classical breakthrough (e.g. Palmer 1991; Dreybrodt 1990; Dreybrodt et al. 2005; Kaufmann 2002): The initial enlargement is focussed to the very first part of the con- duit, where a typical funnel shape evolves (red solid lines in Fig. 4 top). The calcium concentration increases rap- idly over a distance of just a few meters to attain values around 90 % of the equilibrium value. Most of the con- duit thus only slowly enlarges due to the high-order ki- netic rate law active here. With time, the entrance funnel migrates further into the conduit, and first-order kinetics moves forward. Once the first-order kinetics reaches the exit of the conduit, the conduit starts enlarging almost uniformly over its entire length (blue dashed lines in Fig. 4 top), and flow has become turbulent. The transi- tion from high-order to first-order kinetics at the exit of the conduit characterises this two-fold evolution, and the time first-order kinetics reaches the exit is termed break- through time (TB) in the literature (e.g. Dreybrodt 1996). This breakthrough event occurs at around tB ~ 140,000 years in the case of the water-table conduit in limestone. Gypsum The temporal evolution of a conduit in gypsum is shown in Fig. 5 (top). For the given temperature, the calcium equilibrium concentration is around ceq 15.35 mol/m 3. The evolution of the water-table conduit is essentially the same as in the case of limestone: A preferential enlarge- ment to a funnel shape at the entrance during the early phase (red solid lines in Fig. 5 top), with first-order kinet- ics only active in the first few tens of meters. The high- order kinetics active along the entire remaining part of the conduit ensures a bottleneck for flow, as the conduit width close to the exit remains low and thus flow out is limited. Once the first-order kinetics migrates through the conduit and reaches the exit, the breakthrough event occurs, and from then on the conduit enlarges at con- stant pace (blue dashed lines in Fig. 5 top) under tur- bulent flow conditions. Note, however, that the break- through time for the water-table conduit in gypsum is with tB~24,000 years an order of magnitude smaller than in the case of limestone! Anhydrite The temporal evolution of a conduit in anhydrite is shown in Fig. 6 (top). For the given temperature, the calcium GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV ACTA CARSOLOGICA 46/2–3 – 2017 207 equilibrium concentration is around ceq 45 mol/m 3. At a first glance, the evolution of the conduit width is not too different from the two previous cases. Before break- through, a funnel shape at the entrance part evolves, and the remainder of the conduit widens only very slowly (red solid lines in Fig. 6 top). Once the breakthrough event occurs, the conduit grows at almost uniform pace (blue dashed lines) under turbulent flow conditions. However, two points are very different: (i) The break- through time is tB ~ 500 years, which is very short. (ii) The calcium concentration profiles along the conduit differ from the previous two cases; while calcium con- centration increases rapidly over the first few tens of meters, it remains below 90 % of ceq for the remainder of the conduit for all times before breakthrough. This significant undersaturation is a result of the anhydrite calcium flux rate (see Fig. 3), which is non-linear for the entire dissolution branch. No first-order kinetics is present here, and thus the high-order kinetics keep the undersaturation with respect to calcium over the entire conduit length, and is therefore responsible for the fast evolution. DISSOLUTION IN DEEP SETTINGS In this second part of the results section, we discuss the deep-phreatic conduit, which extends to considerable depth zmax. This depth-extent of the conduit causes the solution to be subject to elevated pressure and tempera- ture conditions, thus changes in the chemistry of the so- lution. fig. 4: Evolution of single lime- stone conduit with time. Shown are conduit width and calcium concentration for several time steps (see notation on lines in years). Red lines indicate period before breakthrough, blue lines after breakthrough, and the flow condition is either laminar (solid lines) or turbulent (dashed lines). top: Limestone water-table con- duit. Bottom: Limestone deep- phreatic conduit. FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES ACTA CARSOLOGICA 46/2–3 – 2017208 Limestone This deep-phreatic conduit (Fig. 4 bottom) evolves in a very different manner than the water-table con- duit. During the early evolution the entrance part is enlarged to a funnel shape as before, but with increas- ing depth the conduit evolution becomes inhibited. The reason for the slow enlargement in greater depth is the decrease in calcium equilibrium concentration: At 50 m depth, the temperature is t(50 m)=11.25 °C, hydraulic pressure is pw(50 m)=0.5 MPa, and thus the calcium equilibrium concentration reduces to ceq(50 m) 2.05 mol/m 3, which is around 3 % lower than the surface value. This lower calcium equilibrium con- centration results in slower enlargement with depth, while the ascending part of the conduit (last half of the conduit) enlarges more, as here the calcium equi- librium concentration increases again and allows for faster dissolution (red solid lines). With time, first- order kinetics will also be established in this case, and from then on a breakthrough event occurs, and after that the conduit grows with almost uniform enlarge- ment rates. Note that we increased the initial conduit width to di=1.5 mm for the deep-phreatic conduit to achieve a breakthrough time comparable to the water- table conduit (see Kaufmann et al. 2014, for further details). Note that for a significant depth extent of the deep- phreatic conduit, the elevated temperatures in that depth may result in a reduced calcium equilibrium concentra- tion below the actual calcium concentration, thus pre- fig. 5: As fig. 4, but for gypsum conduit. GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV ACTA CARSOLOGICA 46/2–3 – 2017 209 cipitation would occur and the conduit will clog. We will discuss this clogging process later. Gypsum For the deep-phreatic conduit in gypsum we again have increased the initial conduit width to dini=1.6 mm to ob- tain a similar breakthrough time (tB ~ 30, 000 years). The evolution of this deep-phreatic conduit is, however, very different from the evolution of a similar conduit in lime- stone (Fig. 5 bottom): Before breakthrough, we can dis- tinguish three different regions. (i) An area close to the entrance of the conduit (100−200 m into the fracture), where the first-order kinetics results in a funnel-shape enlargement of the entrance part. (ii) An area of rela- tively uniform growth (200−1500 m into the conduit), which has already enlarged to the centimeter-scale be- fore breakthrough. (iii) The area around the exit of the conduit, which is still small, enlarging at a very slow pace and thus resulting in the bottleneck for flow responsible for keeping the calcium concentration in the high-order regime over large parts of the conduit. The reason for the different evolution before breakthrough is the pro- grade dependence of the calcium equilibrium concentra- tion for gypsum for temperatures below 30 °C: At 50 m depth, the calcium equilibrium concentration increases to ceq(30 m) 15.61 mol/m 3, which is around 2 % larger than the surface value. Thus while the majority of the conduit experiences high-order kinetics before break- fig. 6: As fig. 4, but for anhydrite conduit. FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES ACTA CARSOLOGICA 46/2–3 – 2017210 through, in the lower parts of the conduit the solution is slightly more aggressive due to the increase in ceq. This is responsible for the widespread enlargement in the deeper parts of the deep-phreatic conduit. After breakthrough is established, the entire con- duit enlarges at almost constant pace (blue dashed lines in Fig. 5 bottom) under turbulent flow conditions. Anhydrite The above mentioned significant undersaturation with respect to calcium also changes the evolution for an anhy- drite conduit reaching deeper into the subsurface. Here, no real difference from the water-table conduit evolu- tion is observed, and breakthrough times for both flow paths are very similar. The drop in calcium equilibrium concentration due to the elevated temperature and water pressure at 30 m depth resulting in a calcium equilibrium concentration of ceq(30 m) ~ 44.08 mol/m 3 and thus a re- duction of around 2 % relative to the surface value has no significant effect on the strong undersaturation in the conduit, the high-order kinetics of the calcium flux rate for anhydrite dominates the evolution by far. We note, however, that at this point we neglected the precipitation of gypsum, which will occur in the an- hydrite conduit due to the differences in calcium equi- librium concentration, which will of course change the conduit evolution. We come back to this point in the next section. fig. 7: Clogging of single conduit with time by precipitation. Note that only the first 200 m of the conduits is shown. Evolution of width and calcium concentration for limestone conduit (top) and for gypsum conduit (bottom). GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV ACTA CARSOLOGICA 46/2–3 – 2017 211 DISSOLUTION AND PRECIPITATION IN SHALLOW SETTINGS In this third part of the results section, we will look into the problem of clogging the conduit by precipitation of the corresponding mineral. We will focus first on lime- stone and gypsum as soluble rocks, and treat anhydrite separately because of the precipitation of gypsum in an anhydrite conduit. We argue that the conduit is now subject to inflow with supersaturated solution with respect to calcium. An example would be water passing through a gypsum layer, being saturated with respect to calcium, then flowing into a limestone conduit. We restructure the conduit by addressing a length of l = 2000 m and an initial width of dini = 2.0 mm, and applying the hydraulic head difference of ∆h = 10 m. Limestone We first consider limestone as conduit material, and provide a solution to the conduit, which is slightly su- persaturated with respect to calcium: cin= 2.13 mol/m 3, ceq = 2.12 mol/m 3. In Fig. 7, the evolution of this con- duit is shown. The conduit width starts reducing along the entrance section of the conduit, as here the (small) super-saturation is largest. The deposition becomes smaller along the conduit, as the excess calcium in the solution is consumed. An inverse funnel shape ap- pears, as more and more calcite is deposited, and after tC~856 years the entrance part of the conduit is closed (tC-clogging time) and flow through the conduit stops. Note that still the majority of the initial void volume is present, as the conduit has only been sealed off by a plug. Gypsum We now make up a conduit of gypsum, and feed a solu- tion with cin=15.37 mol/m 3, thus slightly supersaturated when compared to the calcium equilibrium concentra- tion of ceq=15.35 mol/m 3. This conduit closes essentially in the same manner as the limestone conduit, but just in a fraction of the time (TC~9 years). Anhydrite In section 4.1, we have seen that an anhydrite conduit evolves very fast due to the non-linear kinetics of the anhydrite dissolution. However, once the calcium con- centration in the anhydrite conduit reaches the calcium equilibrium concentration of gypsum (~15 mol/m3), gypsum starts to precipitate: We have seen in section 4.2 that a gypsum conduit can become clogged by precipitation of gypsum in only a couple of years, even for small super-saturation with respect to gypsum. In Fig. 8, the change in conduit width (wall retreat) is shown for a system at T=10 °C and both anhydrite and gypsum. The calcium concentration range chosen fo- cuses around the calcium equilibrium concentration for gypsum. While the enlargement of a conduit in anhydrite with a chemical composition around c~15 mol/m3 is with ∆r 6−8 mm/yr rather constant, the situation for pre- cipitation of gypsum dramatically changes the evolution. Once the calcium equilibrium concentration for gypsum ( = 15.3 mol/m 3) is passed, gypsum will start to pre- cipitate. If the calcium concentration in the conduit ex- ceeds just slightly, the deposition rate of gypsum quickly reaches values larger than ∆r>−10 mm/yr, thus outpacing the removal of anhydrite by far. The conduit will clog very soon along the entrance part, and flow is inhibited. The reason for this quick clogging of course is the difference in atomic mass and density for anhydrite and gypsum, which translates into these different retreat rates according to (4). As we have stated in section 3, the precipita- tion flux rates for both anhydrite and gypsum are not well known. We therefore plotted the deposition rate for gypsum also for non-linear flux rate laws (n1=1.5, dashed red line; n1=2, dotted red line). We observe a fairly substantial impact of the non-linearity to gypsum precipitation, which, however, will still outpace the dis- solution of anhydrite. fig. 8: Wall retreat ∆r as a function of calcium concentration c for anhydrite (blue solid line) and gypsum (red lines). for gyp- sum, different rate-equation exponents for precipitation are shown (n1 = 1.0, solid, n1 = 1.5, dashed, n1 = 2.0, dotted). The dashed grey line marks the calcium equilibrium concentration for gypsum at t = 10 °C. FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES ACTA CARSOLOGICA 46/2–3 – 2017212 We speculate that this behaviour will have an im- pact on problems as during the drilling of the hydrother- mal drill hole in Staufen (Sass & Burbaum 2010), which connected two aquifers and caused flow of under-satu- rated solution through an anhydrite lens. As the anhy- drite gets dissolved, the calcium concentration reaches the threshold, when gypsum starts to precipitate. If the gypsum precipitate is able to get firmly attached to the fracture wall, clogging of the anhydrite fracture can oc- cur. However, in the case of very high flow-rate condi- tions, the gypsum precipitate might be flushed out, be- fore it is attached to the fracture wall, thus the fracture remains open. Multi-material conduit In this last section, we pick up the discussion on incep- tion horizons from the introduction. There are numerous observations of faults and bedding planes in limestone more favorable to karstification than others (e.g. Filip- poni et al. 2009; Plan et al. 2009, and references therein). Often, these inception horizons have been covered with pyrite (FeS2), which then oxidized according to (Ritsema & Groenenberg 1993): The first reaction describes the oxidation of pyrite. The latter reaction occurs at the boundary of the fracture wall, where the sulfate reacts with calcite to form gyp- sum. This gypsum, which covers the fracture, can then readily dissolve and provides a rapid initial enlargement of the limestone fracture. Note that the carbon dioxide released can then be dissolved in water and thus increase the calcium equilibrium concentration, water becomes again undersaturated with respect to calcite and can dis- solve additional limestone. We show such an example in Fig. 9. Here, a lime- stone conduit with a length of l = 2000 m and an initial width of dini = 0.5 mm is covered with a gypsum layer of 5 mm thickness. Flow is driven from left to right by the hydraulic head difference of ∆h = 10 m, the incoming so- lution is aggressive (cin = 0 mol/m 3). The conduit experiences a two-stage evolution. First, the gypsum starts dissolving up to its saturation ( 15 mol/m 3). The entrance enlarges as funnel shape, the remainder of the conduit only slowly enlarges. After 10,000 years, the gypsum along the entrance sec- tion of the conduit is gone, the conduit here starts evolv- ing by dissolving limestone with its lower saturation ( 2 mol/m 3). As the solution saturated with respect to limestone is still aggressive with respect to the gyp- sum, the fast gypsum dissolution is pushed further into the conduit. This scenario is similar to the limited dissolution proposed by Romanov et al. (2010) and Gabrovšek and Stepišnik (2011). If there is limitation of soluble mate- rial in the direction perpendicular to flow, enlargement becomes more uniform along the entire conduit. In our case, after 10,000 years the gypsum vanishes along the entire conduit, flow rate increases and dissolution of the exposed limestone then accelerates. When we compare breakthrough times of this gypsum/limestone conduit (tB~21,000 years) to that of a similar conduit only composed of limestone fig. 9: Evolution of limestone conduit with gypsum coating. GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV ACTA CARSOLOGICA 46/2–3 – 2017 213 (tB~140,000 years, see Fig. 4), we find an acceleration of evolution by an order of magnitude. This faster evolu- tion can explain the importance of an inception horizon, which is chemically distinct from the other conduit by the gypsum cover. The inception horizon evolves way faster than the surrounding conduit and thus captures flow and provides a preferential pathway. The limited dissolution due to the thin thickness of the gypsum cover plays an important role. We have re- duced for the example conduit discussed above the thick- ness of the gypsum cover from 5 to 1 mm, and break- through time becomes even shorter (tB~13,000 years). Thus two mechanisms can be identified as important for preferential karstification along inception horizons: The faster dissolution of the gypsum precipitate, and the more uniform enlargement due to the limited dissolu- tion effect. We can of course also speculate about the opposite effect: A limestone conduit coated with gypsum can clog, if the gypsum layer has a non-uniform thickness. If a thinner gypsum coating within the conduit is removed first and the limestone exposed in that part, the solution arriving is significantly supersaturated with respect to limestone, and limestone will precipitate and can eventu- ally clog flow of the entire conduit. We have developed a formal framework for describing the evolution of a single isolated conduit in a soluble rock, which can be enlarged by dissolution of the rock, and which can clog as a result of precipitation of a min- eral in solution. The single conduit can consist of lime- stone, gypsum, or anhydrite, or a combination of these soluble rock types. As stated in Introduction, we defined three key points for our analysis, which we will answer now: i Shallow and deep flow and evolution in con- duits for different soluble rock types (lime- stone, gypsum, anhydrite): A first result, already established in the literature, is the strong dependence of the evolution time on the type of soluble rock considered: While the enlargement of conduits in limestone under natural hydraulic condi- tions can be between 1000 and 10000 years, conduits in anhydrite and gypsum evolve much faster, on timescales of 10–100 years. Deeper flow paths as considered in the bathyphre- atic evolution of caves depend on the change of environ- mental parameters with depth: Temperature, water and rock pressure generally increase with depth, changing hydraulic properties (e.g. water viscosity), but also the chemical properties (e.g. calcium equilibrium concen- tration). Here, the retro-grade dependence of calcium equilibrium concentration on temperature for limestone and anhydrite are responsible for a reduction of enlarge- ment with depth, thus bathyphreatic flow paths evolve slower, when compared to water-table flow paths. For gypsum, however, the pro-grade dependence of calcium equilibrium concentration on temperature (at least for temperatures below 30 °C) does not inhibit evolution with depth. ii Precipitation in conduits: Once solution becomes supersaturated with respect to calcium, either through changes in temperature and/ or water pressure or evaporation, the conduit can clog due to precipitation. As the precipitation of soluble rock is largest along the inflow part of the fracture, clogging creates a plug inhibiting flow, but keeps large parts of the conduit further downstream still open. iii Evolution of conduits composed of several sol- uble rock layers: If a conduit consists of more than one material, the evolution becomes more complicated. In our example of a limestone conduit coated with gypsum (e.g. from the conversion of pyrite into gypsum) evolves more quickly, when compared to a pure limestone conduit. Here, the gypsum coating is quickly removed along the entire con- duit, thus the remaining limestone part has a larger di- ameter, which for the ongoing evolution exerts a strong control on the time of evolution. CONCLUSIONS FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES ACTA CARSOLOGICA 46/2–3 – 2017214 We would like to thank two anonymous referees, whose comments improved the manuscript substantially. GK acknowledges funding from the DFG under research grant KA1723/6- 2. This project has been carried out at the Karst Research Institute of Slovenia (ZRC SAZU) during a sabbatical stay of GK. Figures were prepared us- ing GMT software (Wessel & Smith 1998). ACKNOWLEDGEMENTS REFERENCES Beek, W. J. & K. M. K. Muttzall, 1975: transport Phenomena.-John Wiley, London. Blount, C. & F. Dickso, 1973: Gypsum-anhydrite equilib- ria in systems CaSO4-H2O and CaSO4-NaCl-H2O.- Am. Mineralogist, 58, 323–331. Bretz, J. H., 1942: Vadose and phreatic features of lime- stone caves.-. J. Geology 50, 675–811. DOI: https:// doi.org/10.1086/625074. Buhmann, D. & W. Dreybrodt, 1985a: The kinetics of calcite dissolution and precipitation in geologi- cally relevant situations of karst areas. 1. Open sys- tem.- Chem. Geol., 48, 189–211. DOI: https://doi. org/10.1016/0009-2541(85)90046-4. Buhmann, D. & W. Dreybrodt, 1985b: The kinetics of calcite dissolution and precipitation in geological- ly relevant situations of karst areas. 2. Closed sys- tem.- Chem. Geol., 53, 109–124. DOI: https://doi. org/10.1016/0009-2541(85)90024-5. Davis, W. M., 1930: Origin of limestone caverns.- Geol. Soc. Am. Bull., 41, 475–628. DOI: https://doi.org/ 10.1130/GSAB-43-663. Dreybrodt, W., 1988: Processes in karst Systems. Springer, 278 pp., Berlin. Dreybrodt, W., 1990: The role of dissolution kinetics in the development of karst aquifers in limestone: a model simulation of karst evolution.- J. Geol., 98, 5, 639–655. DOI: https://doi.org/10.1086/629431. Dreybrodt, W., 1996: Principles of early development of karst conduits under natural and man-made condi- tions revealed by mathematical analysis of numeri- cal models.- Water Resour. Res., 32, 9, 2923–2935, https://doi.org /10.1029/96WR01332. Dreybrodt, W. & D. Buhmann, 1991. A mass transfer model for dissolution and precipitation of calcite from solutions in turbulent motion.- Chem. Geol., 90, 107–122. DOI: https://doi.org/10.1016/0009- 2541(91)90037-R. Dreybrodt, W., Eisenlohr, L., Madry, B. & S. Ringer, 1997: Precipitation kinetics of calcite in the system CaCO3-H2O-CO2: The conversion to CO2 by the slow process H++HCO−3 →CO2+H2O as a rate lim- iting step.- Geochem. Cosmochem. Acta, 61, 18, 3897–3904. DOI: https://doi.org/10.1016/S0016- 7037 (97)00201-9. Dreybrodt, W., Gabrovšek, F. & D. Romanov, 2005: Processes of Speleogenesis: A modeling approach. Založba ZRC, 375 pp., Ljubljana. Dreybrodt, W. & G. Kaufmann, 2007: Physics and chem- istry of dissolution on subaerially exposed soluble rocks by flowing water films.- Acta Carsologica, 36, 3, 357–367. DOI: https://dx.doi.org/10.3986/ ac.v36i3.169. Duan, Z. & D. Li, 2008. An improved model for the calculation of CO2 solubility in aqueous solutions containing Na+, K+, Ca2+, Mg2+, Cl−, and SO2−.- Geo- chem. Cosmochem. Acta, 72, 5128–5145. DOI: https://doi.org/10.1016/j.marchem.2005.09.001. Eisenlohr, L., Meteva, K., Gabrovšek, F. & W. Dreybrodt, 1999: The inhibiting action of intrinsic impurities in natural calcium carbonate minerals to their dissolu- tion kinetics in aqueous H2O-CO2 solutions.- Geo- chim. Cosmochem. Acta, 63, 6, 989–1001. DOI: https://dx.doi.org/10.1016/S0016-7037(98)00301-9. Filipponi, M. & P.-y. Jeannin, 2006: Is it possible to pre- dict karstified horizons in tunneling?- Austrian Journal of Earth Sciences, 99, 24–30. Filipponi, M., Jeannin, P. & L. Tacher, 2009: Evidence of inception horizons in karst conduit networks.- Geomorphology, 106, 86–99. DOI: https://doi.org/ 10.1016/ j.geomorph.2008.09.010. Ford, D. C. & R.O. Ewers, 1978: The development of limestone cave systems in the dimensions of length and depth.- Can. J. Earth Sci., 15, 1783–1798. DOI: http://dx.doi.org/10.5038/1827-806X.10.3.1. https:// doi.org/10.1139/e78-186. GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV ACTA CARSOLOGICA 46/2–3 – 2017 215 Gabrovšek, F. & U. Stepišnik, 2011: On the formation of collapse dolines: A modeling perspective. Geomor- phology, 134, 23–31. DOI: https://doi.org/10.1016/j. geomorph.2011.06.007. Gobran, G. & S. Miyamoto, 1985: Dissolution rates of gypsum in aqueous salt solutions.- Soil Science, 140, 2, 89–93. DOI: https://10.1097/00010694-198- 508000-00002. Groves, C. G. & A. D. Howard, 1994a. Minimum hydro- chemical conditions allowing limestone cave devel- opment. Water Resour. Res. 30 (3), 607–615, https:// dx.doi.org/10.1029/93WR02945. Groves, C. G. & A.D. Howard, 1994b. Early develop- ment of karst systems 1. Preferential flow path enlargement under laminar flow. Water Resour. Res., 30, 10, 2837–2846. DOI: https://dx.doi. org/10.1029/94WR01964. Howard, A. D. & C.G. Groves, 1995. Early develop- ment of karst systems 2. Turbulent flow. Water Resour. Res., 31, 1, 19–26. DOI: https://dx.doi. org/10.1029/94WR01303 James, A. & A.Lupton, 1978. Gypsum and anhydrite in foundations of hydraulic structures. Geotechnique, 28, 3, 249–272. DOI: https://doi.org/10.1680/ geot.1978.28.3.249. Jeppson, R. W., 1976. Analysis of flow in pipe networks. Ann Arbor Sci. Pub., 164 pp., Ann Arbor. Jeschke, A., Vosbeck, K. & W. Dreybrod, 2001. Surface controlled dissolution rates of gypsum in aque- ous solutions exhibit nonlinear dissolution rates. Geochem. Cosmochem. Acta, 65, 1, 27–34. DOI: https://doi.org/10.1016/S0016-7037(00)00510-X Jeschke, A. A., 2002; Lösungskinetik von gips und An- hydrit. Ph.D. thesis, Universität Bremen, Bremen, Germany. Jones, F. & R. Schoonove, 2002. handbook of Mass Mea- surements. CRC Book, 336 pp. Kaufmann, G., 2002: karst conduit evolution. In: Gabrovšek, F. (Ed.), Evolution of karst: from Preka- rst to Cessation. Carsologica, ZRC SAZU, Postojna- Ljubljana, pp. 327–338. Kaufmann, G., Gabrovsek, F. & D. Romanov, 2014: Deep conduit flow in karst aquifers revisited.- Water Re- sour. Res, 50, 6, 4821–4836. DOI: https://dx.doi. org/10.1002/2014WR015314. Kaufmann, G., Romanov, D. & T. Hiller, 2010: Model- ling three-dimensional karst aquifer evolution us- ing different matrix-flow components.- J. Hydrol., 388, 241–250. DOI: https://dx.doi.org/10.1016/j. jhydrol.2010.05.001. Kestin, J., Sokolov, M. & W. Wakeham, 1978: Viscosity of liquid water in the range −8 °C to 150 °C. J. Phys. Chem. Ref. Data, 7, 3, 941–948. DOI: https://dx.doi. org/10.1063/1.555581. Lebedev, A. & A. Lekhov, 1990: Dissolution kinetics of natural gypsum in water at 5–25 °C.- Geochem. Int., 27, 85–94. Li, J. & Z. Duan, 2011: A thermodynamic model for the prediction of phase equilibria and speciation in the H2O-CO2-NaCl-CaCO3-CaSO4 system from 0 to 250 °C, 1 to 1000 bar with NaCl concentrations up to halite saturation.-Geochem. Cosmochem. Acta, 75, 4351–4376. DOI: https://doi.org/10.1016/j.gca. 2011.05.019 Lowe, D., 1992: The origin of limestone caverns: an incep- tion horizon hypothesis. Ph.D. thesis, Manchester Polytechnic, United Kingdom. Lowe, D., 2000: Role of stratigraphic elements in spe- leogenesis: the speleo inception concept. In: Klim- chouk, A., Ford, D., Palmer, A. & W. Dreybrodt (Eds.), Speleogenesis, evolution of karst aquifers. Na- tional Speleological Society, Huntsville (Alabama), pp. 65–76. Lubitz, C., Motagh, M., Wetzel, H.-U. & H. Kaufmann, 2013: Remarkable urban uplift in Staufen im Breis- gau, Germany: Observations from terrasar-x insar and leveling from 2008 to 2011.- Remote Sensing, 5, 3082–3100. DOI: https://dx.doi.org/10.3390/rs- 5063082. Mehrbach, C., Culberson, C.H.and Hawley, J. & R. Pytko- wicz, 1973: Measurement of the apparent dissocia- tion constants of carbonic acid in seawater at atmo- spheric pressure.- Limnology and Oceanography, 18, 6, 897–907. DOI: https://dx.doi.org/10.4319/ lo.1973.18.6.0897. Millero, F. J., 1979: The thermodynamics of the car- bonate system in seawater.- Geochemica et Cos- mochemica Acta 43, 1651–1661, https://dx. doi. org/10.1016/0016-7037(79)90184-4. Mucci, A., 1983: The solubility of calcite and aragonite in seawater at various salinities, temperatures, and one atmosphere total pressure. -American Jour- nal of Science, 283, 781–799. DOI: https://dx.doi. org/10.2475/ajs.283.7.780. Palmer, A. N., 1991: Origin and morphology of lime- stone caves.- Geol. Soc. Am. Bull., 103, 1–21. DOI: https://doi.org/10.1130/0016-7606(1991)103%3C0 001:OAMOLC%3E2.3.CO;2. FRACTURE EVOLUTION IN SOLUBLE ROCKS: FROM SINGLE-MATERIAL FRACTURES TOWARDS MULTI-MATERIAL FRACTURES ACTA CARSOLOGICA 46/2–3 – 2017216 Parkhurst, D. & C. Appelo, 2013: Description of input and examples for PHREEqC version 3 – A com- puter program for speciation, batch-reaction, one- dimensional transport, and inverse geochemical calculations.- U.S. Geological Survey Techniques and Methods, Book 6, Chap. A43, 497 p., available online http://pubs.usgs.gov/tm/06/a43 (Accessed October 10th, 2017) Plan, L., Filipponi, M., Behm, M., Seebacher, R. & P. Jeutter, 2009: Constraints on alpine speleogenesis from cave morphology – a case study from the east- ern Totes Gebirge (Northern Calcareous Alps, Aus- tria).- Geomorphology, 106, 118–129. DOI: https:// doi.org/10.1016/j.geomorph.2008.09.011. Plummer, L. N., Wigley, T. M. L. & D.L. Parkhurst, 1978: The kinetics of calcite dissolution in CO2–water sys- tems at 5 °C to 60 °C and 0.0 to 1.0 atm CO2.- Am. J. Sci., 278, 179–216. DOI: htpps://doi.org/ 10.2475/ ajs.278.2.179 Rhoades, R. & N.M. Sinacori, 1941; Patterns of ground- water flow and solution.- J. Geology, 49, 785–794. DOI: https://doi.org/10.1086/625014. Ritsema, C. & J. Groenenberg, 1993: Pyrite oxida- tion, carbonate weathering, and gypsum forma- tion in a drained potential acid sulfate soil.- Soil Science Society of America Journal, 57, 968–976. DOI:https://dl.sciencesocieties.org/publications/ss- saj/abstracts/57/4/SS0570040968. Romanov, D., Kaufmann, G. & T. Hiller, 2010: Karstifica- tion of aquifers interspersed with non-soluble rocks: from basic principles towards case studies.- Eng. Geol., 116, 261–273. DOI: https://doi.org/10.1016/j. enggeo.2010.09.008. Sass, I. & U. Burbaum, 2010: Damage caused to the his- toric town of Staufen (Germany) caused by geo- thermal drillings through anhydrite-bearing for- mations.- Acta Carsologica, 39, 2, 233–245. DOI: http://dx.doi.org/10.3986/ac.v39i2.96. Svensson, U. & W. Dreybrodt, 1992: Dissolution kinet- ics of natural calcite minerals in CO2-water systems approaching calcite equilibrium. Chem. Geol., 100, 129–145. DOI: https://doi.org/10.1016/0009-2541- (92)90106-F. Swinnerton, A. C., 1932: The origin of limestone cav- erns.- Geol. Soc. Am. Bull., 34, 662–693. DOI: https://doi.org/10.1130/GSAB-43-663. Szymczak, P. & A. Ladd, 2006: A network model of chan- nel competition in fracture dissolution.- Geophysi- cal Research Letters, 33, L05401. DOI: https://doi. org/10.1029/2005GL025334. Szymczak, P. & A. Ladd, 2009. Wormhole formation in dissolving fractures.- J. Geophys. Res. M 114, B06203. DOI: https://doi.org/10.1029/2008JB006122. Szymczak, P. & A. Ladd, 2011: The initial stages of cave formation: Beyond the one-dimensional paradigm. Earth and Planetary Science Letters, 201, 424−432. DOI: https://doi.org/10.1016/j.epsl.2010.10.026. Weiss, R. F., 1974: Carbon dioxide in water and seawater: the solubility of a non-ideal gas.- Marine Chemis- try 2, 203–215, https://doi.org/10.1016/0304-4203- (74)90015-2. Wessel, P. & W. H. F. Smith, 1998. New, improved ver- sion of generic mapping tools released.- EOS, 79, 579. DOI: https://doi.org/10.1029/98EO00426 Wissbrun, K., French, D. & A. Patterson Jr., 1954: The true ionization constant of carbonic acid in aqueous solution from 5 to 45°.- J. Phys. Chem., 58, 9, 693– 695. DOI: htpps://doi.org/10.1021/j150519a004. Worthington, S. R. H., 2001: Depth of conduit flow in unconfined carbonate aquifers.- Geology, 29, 4, 335–338. DOI: https://doi.org/10.1130/0091-7613- (2001)029%3C0335:DOCFIU%3E2.0.CO;2. GEORG KAUFMANN, FRANCI GABROVŠEK & DOUCHKO ROMANOV