Carbonates Evaporites DOI 10.1007/s13146-017-0337-6

ORIGINAL ARTICLE

Mathematical modeling of karstogenesis: an approach based on fracturing and hydrogeological processes Antoine Lafare1 • Herve´ Jourde2 • Ve´ronique Le´onardi2 • Se´verin Pistre2 Nathalie Do¨rfliger3

•

Accepted: 29 January 2017 Ó Springer-Verlag GmbH Germany 2017

Abstract The karstogenesis of a synthetic aquifer is analyzed as a function of the matrix permeability and the fracture density affecting the carbonate reservoir. The synthetic carbonate aquifer generation results from numerical simulations based on Discrete Fracture Network (DFN) and groundwater flow simulations. The aim is to simulate karstogenesis processes in an aquifer characterized by a fracture network, while matching field reality and respecting geometrical properties as closely as possible. DFN are simulated with the soft REZO3D that allows the generation of 3-D realistic fracture networks, especially regarding the relative position of fractures that control the overall network connectivity. These generated DFN are then meshed and considered to perform groundwater flow simulation with the model GroundWater (GW). At each time step and for each fracture element, flow velocity and the mean groundwater age are extracted and used in an analogical dissolution equation that computes the aperture evolution. This equation relies on empirical parameters calibrated with former speleogenesis studies. In this paper, karstogenesis simulations are performed using fixed-head hydraulic boundary conditions within a single stratum as a function of two fracture density settings. The results are

& Antoine Lafare [email protected] 1

British Geological Survey, Nicker Hill, Keyworth, Nottingham NG12 5GG, UK

2

Laboratoire Hydrosciences Montpellier, UMR 5569 CNRSIRD-UM1-UM2, Universite´ Montpellier 2, Place Euge`ne Bataillon, 34095 Montpellier, France

3

Water, Environment and Ecotechnologies Division, BRGM, 3 Avenue Claude-Guillemin, BP 36009, 45060 Orle´ans Cedex 02, France

interpreted in terms of head fields, mean groundwater age distributions; while total flow rate and average aperture are analyzed as a function of time. The effect of fracture density and rock matrix permeability are then assessed. Keywords Speleogenesis  Karst  Fracturing  Hydrogeology  Modeling  Groundwater

Introduction Karst hydrogeological systems are difficult to manage because of their great heterogeneity and vulnerability to contaminations. Indeed, such aquifers exhibit a dual flow system (Kiraly et al. 1998) with a highly transmissive conduits network, and a capacitive bulk mass of fractured rock. With the goal of managing karst groundwater resource either for quantity or quality issues, it is necessary to use adequate numerical models to simulate various scenari. Models used so far, such as lumped, reservoirs or distributed models (Mangin 1975; Sauter 1992; Do¨rfliger et al. 2009; Fleury et al. 2009; Tritz et al. 2011; Mazzilli et al. 2012), either do not take into account the geometry of the karst drainage system or account for a simplified geometry. A better knowledge of this geometry could improve the management of such aquifers, specifically for quality or vulnerability issues and mass transfer modeling. Several numerical modeling approaches were undertaken to simulate the processes of karst conduit genesis (Palmer 1991; Dreybrodt 1996; Kaufmann and Braun 1999, 2000; Bauer 2003; Liedl et al. 2003; Dreybrodt et al. 2005). These existing methods are mainly based on the physical and chemical laws driving the carbonate dissolution process (taking into account calcite saturation of the water and partial pressure of carbon dioxide). As a

123

Carbonates Evaporites

consequence, these previous works provide a well-documented knowledge on the kinetics of the carbonate dissolution processes in karst systems. Nevertheless, these models are mainly applied on simplified initial void networks, which do not match the fracturing and geological reality. However, the initial geometry of the void network (namely the geological structures such as fractures and bedding planes, and their connectivity) determines an important part of the final pattern of the speleological network (Klimchouk and Ford 2000; Bakalowicz 2006; Ford and Williams 2007). Considering realistic simulated fracture sets as initial void networks could therefore improve the understanding of speleogenetic processes. In this aim, a numerical model is developed (Lafare et al. 2009; Lafare 2011), which involves a pseudo-statistic fracturing generator [REZO3D, (Jourde 1999; Josnin et al. 2002; Jourde et al. 2002, 2007)], coupled to a finite element groundwater simulator [GroundWater, (Cornaton and Perrochet 2006a, b; Cornaton 2007; Cornaton et al. 2008)]. The modeling of the genesis of the karst drainage system is based on an analogical empirical polynomial equation considering the pore velocity and mean age of water as main parameters. The polynomial form of this equation was designed taking the mathematical bases of the MOPOD model (Bloomfield and Barker 2005; Bloomfield et al. 2005) as reference. The simulations are performed on the basis of a time step, whose duration depends on the simulated scenario. The mean age of the water is used to represent the progressive decrease of the chemical dissolving potential of the water in contact with the carbonate rock within the aquifer (Dreybrodt and Eisenlohr 2000). The first simulator—REZO3D—allows generating three-dimensional networks constituted by plane fractures, whose lengths and spatial distribution respect mechanical and statistical laws. These networks are then processed to write finite element meshes which constitute the basis of further groundwater flow and transport simulations, considered in the karstogenesis processes. The polynomial coefficients of the analogical empirical equation for the karstogenesis were adjusted to approach the results (evolution of fracture aperture pattern and of the total flow rate with time) obtained by the reference simulations based on the dissolution kinetic laws of carbonate (Palmer 1991; Dreybrodt 1996; Palmer 2000; Dreybrodt et al. 2005). Once these calibration procedures were carried out, the methodology was applied to more complex situations (Lafare 2011). In this paper, simulations involving fracture networks generated with REZO3D are presented. These networks consist in two orthogonal families of 2D rectangular fractures embedded in a single stratum represented by a 3D rocky matrix. Karstogenesis is simulated and analyzed according to two settings characterized by different fractures densities. The influence of the fracture density on

123

the karstogenesis processes is therefore explored, in terms of evolution of the distribution of the groundwater age, the head field, and the total flow rate. Besides, the effect of the hydrodynamic characteristics of the carbonate matrix on the synthetic karst network development is also assessed, using three simulation scenarios characterized by various matrix hydraulic conductivities.

Materials and methods The main steps of the methodology are summarized in the diagram presented in Fig. 1, and extensively described in Lafare et al. (2009) and Lafare (2011). A first set of preprocessing procedures (Fig. 1, dotted lines) aims to generate an initial void network using the fracture generator REZO3D, and to process it to create a finite element mesh suitable for groundwater flow and transport computations. Statistical and mechanical parameters controlling the generation and distribution of the fracture networks are selected according to the characteristics of the fracture network that is to be modeled (fracture density, length, number of layers, etc.). Two orthogonal fracture families are then produced, each of these families being characterized by the mean and standard deviation of the fractures length, the density represented by the number of initiation points, the number (one in the presented cases), and average thickness of sedimentary layers. Two sets of parameters also control the mechanical relationship between fractures, either of the same family or not. These parameters determine the shape of a relaxation zone around the fractures, preventing the initiation or the development of further fracture (Jourde 1999; Josnin et al. 2002). The meshing procedure then allows the generation of a fractured reservoir system composed of different flow elements. According to the meshing option chosen, 4-noded fracture elements and/or 2-noded pipe elements are embedded in a comparatively low hydraulic conductivity matrix made of homogeneous 8-noded brick elements. The initial apertures of each fracture element, as well as the hydrodynamic parameters of the matrix, are then set by a dedicated procedure. The main simulation procedure (Fig. 1, bold lines) consists in a computing loop involving flow and mean age computations with GroundWater (GW) and a dissolution procedure that modifies the aperture of each fracture and/or conduit element according to the results of the previous computations. For the groundwater flow simulations, a matrix continuum is considered and a discrete fracture approach is used. The Hagen–Poiseuille law controls the 2D fracture elements hydraulic conductivity under saturated and laminar flow conditions, while under turbulent flow conditions the Manning–Strickler law is used. GW is

Carbonates Evaporites Fig. 1 Outline of the main steps of the modeling approach

also able to compute the mean groundwater age within the synthetic aquifer, defined as a relative quantity with respect to a starting location where age is supposed to be zero (Cornaton 2003; Cornaton and Perrochet 2006b; Cornaton 2007). For each time step, at the end of the GW execution, the distributed flow velocities and groundwater ages are extracted for each fracture element. Then, a ‘‘dissolution’’ procedure performs a loop on each element, incorporating the previously extracted variables in the aforementioned analogical empirical polynomial dissolution equation, to compute a new aperture. The next simulation time step can then be launched using the new apertures. The analogical equation (Eq. 1) was elaborated considering polynomial functions of the magnitude of the flow velocity, and of the mean groundwater age: dai avbi ¼ dt 1 þ cAdi

i ¼ 1; . . .; N

ð1Þ

where a, b, c, and d refer to the polynomial coefficients allowing the adjustment of the equation, ai [L], vi [LT-1] and Ai [T] represent, respectively, the aperture, the magnitude of the flow velocity and the mean groundwater age for the element i. The groundwater age is considered here as an enlargement inhibitor factor (as the time spent by the water within the aquifer increases, its potential aggressiveness decreases) and is then placed at the denominator. Indeed, as the time spent by the water within the carbonate aquifer increases, its amount of dissolved carbonate increases while its potential aggressiveness decreases.

Several suitable sets of empirical polynomial coefficients a, b, c, and d were adjusted (Lafare 2011) to approach the results obtained by the reference simulations aiming to simulate the aperture development of a single fracture, based on the calcite dissolution kinetics laws (Palmer 1991; Dreybrodt 1996; Dreybrodt et al. 2005). In this aim, several benchmark simulations on a single one-dimensional fracture were performed, involving simulation design and parameters as close as possible to the works taken as references (Lafare 2011). First, acceptable ranges of polynomial coefficients for each setting were determined using a Monte Carlo approach (random sampling using a uniform distribution) aiming to minimize the error between the simulated dissolution pattern and the reference (Root Mean Square Error, RMSE). Then, sensitivity analyses were performed on these ranges to find the most suitable sets of coefficients. The set that has been chosen for this study is the following: a = 0.0008, b = 0.03, c = 1.0, d = 0.2. This set of coefficients was determined as suitable to deal with planar 2-D fractures as tested by Dreybrodt et al. (2005).

Simulations settings The aim of this study is to explore the influence of two parameters on the karstogenesis of a synthetic karst aquifer: on one hand the influence of the fracture density, and on the other hand the influence of the carbonate matrix permeability.

123

Carbonates Evaporites

The fracture networks generated are based on a same basic geometrical setting consisting in a 100 m long (y direction), 50 m large (x direction), and 2.5 m high single stratum (Fig. 2; Table 1). A first fracture network (Network 1, dense, Fig. 2a) is generated, containing a first family (I) of 200 fractures (on the y direction, with a mean length = 50 m) and a second (II) of 400 fractures (on the x direction, with a mean length = 5 m). The second network (Network 2, Fig. 2b) is generated using the same parameters, except the number of fractures: only 100 fractures for the first family and 200 for the second one. Both networks present orthogonal families of vertical fractures (Fig. 2). The flow in these DFN is driven by a fixed-head boundary condition, with a fixed H = 5 m head at the inlet plane y = -50 m, and a H = 0 m head condition at the outlet plane y = 50 m (see Fig. 2; Table 2). All the fracture elements have a homogeneous initial aperture of 0.01 cm. The simulations are performed for a time period of 10,000 years [a duration that can be sufficient for fissures to grow to traversable size under certain conditions (Palmer 2001)], with 20 time steps lasting 500 years each, and the enlargement equation parameterized as indicated in Table 2. The carbonate matrix is represented by 3-D elements characterized by a hydraulic conductivity Km = 10-6 m/s. To assess the influence of different matrix hydraulic

conductivities, two additional simulations are performed. The matrix hydraulic conductivity rises to Km = 10-5 and then to Km = 10-4 m/s (see Table 2).

Results The results concerning the influence of the fracture density are presented in a first part, with maps representing the evolution of groundwater ages distribution (Figs. 3, 4) and charts showing the temporal evolution of the mean aperture over the length of the fracture element and of the total flow rate (Fig. 5). In a second part, the results for a variable matrix hydraulic conductivity are provided with maps of the evolution of the head fields (Fig. 6) and charts showing the temporal evolution of the mean aperture (Fig. 7) and of the total flow rate (Fig. 8). Influence of the fracture density on the synthetic karst aquifer genesis Temporal evolution of the groundwater age distribution The temporal evolution of the groundwater age distribution is represented for both networks by six snapshots drawn at

Fig. 2 Simulations settings for the fracture networks: Network 1 (a) and Network 2 (b)

Table 1 Values of the parameters used to generate the fracture network 1 and 2

Network 1 (600 fractures) Domain width x (m)

50

50

Domain length y (m)

100

100

Number of layers

1

1

Layer thickness (m)

2.5

2.5

Number of fractures I

200

100

Mean length fractures I (m)

50

50

Number of fractures II

400

200

Mean length fractures II (m)

5

5

I refers to the first fracture family, while II represents the second one

123

Network 2 (300 fractures)

Carbonates Evaporites Table 2 Values of the parameters used within the flow, mean age, and enlargement computations

Network 1 (600 fractures)

Network 2 (300 fractures)

Input hydraulic head (m)

5

5

Output hydraulic head (m)

0

0

Initial aperture a(t=0) (cm)

10-2 -6

Km (m/s)

10

Total duration (years)

10,000

10-2 10-6, 10-5, 10-4 10,000

Number of time steps

20

20

a

0.0008

0.0008

b

0.03

0.03

c

1.0

1.0

d

0.2

0.2

Fig. 3 Groundwater age distribution at six times for the dense network (Network 1)

1000, 2500, 4000, 5500, 7000, and 8500 years (Fig. 3 for Network 1 and Fig. 4 for Network 2, respectively). To emphasize the propagation of relatively ‘‘young’’ water, the age scale presents values from 0 to more than 0.05 years in the legends of Figs. 3 and 4. It appears that the temporal evolution is relatively similar for both fracture networks. At 1000 years, almost the entire aquifer drains water older than 0.05 years except at locations very close to the infiltration limit, where the flow

is not fast enough to allow the young water to penetrate within the system. At 2500 years, young water penetrates a few meters within the fractures directly connected to the input boundary through the initiated preferential enlargement. At 4000 years, the young water penetrates deeper within the aquifer, through the most enlarged and connected flow paths. In the case of the less densely fractured network (Fig. 4), it can be noted that the penetration is more straightforward and uses a lower number of paths

123

Carbonates Evaporites

Fig. 4 Groundwater age distribution at six times for the less dense network (Network 2)

Fig. 5 Temporal evolution of the mean fracture aperture (a) and the total flow rate (b) for both simulations characterized by different densities of fracturing

than in the case of the more densely fractured network (Fig. 3). Accordingly, the young water can then flow slightly farther from the input limit. At 5500 years, the progressive enlargement of the fractures allows the penetration of young (and therefore more aggressive) water down to the discharge limit. The

123

enlarged flow paths converge, or split according to the connectivity of the fracture networks. A characteristic time t, referred to as ‘‘Breakthrough time’’ can be defined as the moment when the ratio between the total flow rate at t and the initial total flow rate reaches some predetermined value, e.g., 103 (Dreybrodt et al. 2005). This increase in

Carbonates Evaporites

Fig. 6 Groundwater hydraulic head distribution at six times for Network 2 with a matrix hydraulic conductivity Km = 10-6 m/s

flow rate is sustained by a flow morphology that takes place from the input to the output limit (Jeannin 1996; Dreybrodt et al. 2005). On the last snapshots (7000 and 8500 years), the enlargement continues and an increasing amount of young water diffuses within the carbonate matrix from the most enlarged fractures. Indeed, the fixed-head boundary condition chosen does not limit neither the flow rate, nor the recharge within the aquifer. As a consequence, steep hydraulic gradients propagate into the stratum towards the discharge limit, between the enlarged fracture elements and the adjacent carbonate matrix. The morphologies obtained (Figs. 3, 4) can be compared to natural ones described for example by Palmer (1991). In the mentioned paper, the main types of cave patterns are described, as well as their relationship to the different types of recharge and porosity that can be encountered in karst systems. In the presented study, the initial porosity is represented by orthogonal fracture networks. As a consequence, ‘‘angular’’ passages are generated. Different morphologies can be isolated according to the different snapshots. Taking the example of the most dense network (Fig. 3), one can see that at 5500 years a limited number of pathways bring the young water to the discharge limit, fed by several branches. This morphology is designated as a

‘‘branchwork’’. On the other hand, at 7000 years the simultaneous enlargement of numerous competing openings is observed, indicating the evolution towards a ‘‘maze’’. According to Palmer (1991), such morphology can be caused by steep hydraulic gradients and low saturation ratios. This is consistent with the fact that at this moment of the evolution, the steep hydraulic gradient propagates within the whole stratum. The role of the recharge is therefore emphasized: if the boundary conditions were different, for example involving a limited flow rate, the results could be different in terms of morphology. Temporal evolution of the mean aperture and the total flow rate The temporal evolutions of the mean aperture of the fractures (Fig. 5a) and of the total flow rate (Fig. 5b) do not seem to depend on fracture density as these evolutions as a function of time are similar for both simulations. Based on the flow rate evolution as a function of time, it can be observed that the breakthrough Time (determined here for Q = 103 9 Q(t=0)) occurs a few years earlier in the case of the less dense network. This is consistent with the previous observation indicating that the penetration of the

123

Carbonates Evaporites

Fig. 7 Groundwater hydraulic head distribution at six times for Network 2 with a matrix hydraulic conductivity Km = 10-4 m/s

a

b

Fig. 8 Temporal evolution of the mean fracture aperture (a) and the total flow rate (b) for three simulations characterized by increasing matrix hydraulic conductivities: 10-6 m/s, 10-5 m/s (med), 10-4 m/s (high)

young water is more direct and uses a lower number of paths for the less densely fractured network. Therefore, the discharge limit can be reached slightly earlier and the flow morphology allowing the breakthrough can settle. This threshold is attained between 5000 and 6000 years, which

123

is consistent with the observation described in the previous ‘‘Temporal evolution of the groundwater age distribution’’: the breakthrough is thought to occur in the time period comprised between the snapshots taken at 5500 years and 7000 years.

Carbonates Evaporites

Influence of the matrix hydraulic conductivity on the evolution of the synthetic karst aquifers Temporal evolution of the groundwater hydraulic head distribution Figures 6 and 7 show the temporal evolution of the groundwater hydraulic head distribution for the simulation scenarios with the lowest (Km = 10-6 m/s) and the highest (Km = 10-4 m/s) matrix hydraulic conductivity, respectively. As for the groundwater age in the ‘‘Temporal evolution of the groundwater age distribution’’, six snapshots are provided from 1000 to 8500 years. It has to be noted that in this section, only the less densely fractured network is considered: the focus is directed on the influence of the carbonate matrix hydraulic conductivity. It can be observed that the shapes of the obtained head fields are very similar, but do not appear at the same time, which indicates a slightly different evolution as a function of time. At 1000 years as well as at 2500 years, the hydraulic gradient is almost uniform on the whole aquifer in the case of the highest matrix hydraulic conductivity (Fig. 7). On the other hand, the hydraulic gradient already appears perturbed at 2500 years in the case of the less permeable matrix (Fig. 6). Such a perturbation doesn’t occur until 4000 years for the more permeable matrix (Fig. 7). At 5500 years, a matured flow morphology can be observed within the aquifer characterized by the less permeable matrix (Fig. 6) while this feature appears only at 7000 years for the more permeable matrix (Fig. 7). It appears that the most permeable rock matrix controls the main part of the flow during a longer time than the less permeable one: the fractures have to attain a higher aperture to gain more control on the whole flow within the aquifer. In the case of a low permeable carbonate matrix, the contrast of permeability between the enlarged fracture and the matrix becomes rapidly considerable. As a consequence, the hydraulic gradient is influenced by the location of the enlarged pathways earlier than for the case of a highly permeable matrix. Temporal evolution of the mean aperture and total flow rate The temporal evolutions of the mean aperture of the fractures (Fig. 8a) and of the total flow rate (Fig. 8b) are then presented. It appears clearly that these evolutions differ according to the simulation setting, even if the shapes are relatively similar. Indeed, the mean aperture of the fractures increases faster carbonate matrix hydraulic conductivity increases (Fig. 8a): a mean aperture of 1 cm is attained at about 6000 years for Km = 10-4 m/s,

6500 years for Km = 10-5 m/s, and 7000 years for Km = 10-6 m/s. The higher hydraulic conductivity of the matrix, combined with a fixed-head boundary condition induces an increased inflow of potentially aggressive water at the beginning of the simulation. In consequence, relatively young water is able to reach a more important number of fracture elements that are then enlarged, increasing the mean fracture aperture. Concerning the evolution of the total flow rate (Fig. 8b), it can be observed that the flow rates are already different at the beginning of the simulations: the fixed-head boundary condition added to an increased whole permeability implies an increased flow rate. Regarding the temporal evolution, during the first times, the flow rate remains almost constant until the enlarged fractured flow paths attain a sufficient aperture to control an appreciable part of the flow (or until the breakthrough occurs). This change occurs later in the case of a more permeable matrix: before 2000 years for Km = 10-5 m/s, 2500 years for Km = 10-4 m/s. The curves then follow a similar slope, the enlarged fracture network controlling an increasing part of the flow. These results in terms of the influence of the carbonate matrix hydraulic conductivities are consistent with the findings issued from studies focusing on the development of a single fracture embedded in a continuous matrix. The breakthrough times are reduced in the case of more permeable matrices, ‘‘the porous continuum acting as a sink’’ (Kaufmann 2003).

Conclusion The karstogenesis modeling is performed for two realistic fractures networks characterized by different fractures density. The realistically irregular geometries of both networks allow obtaining preferential flow paths determined by the connectivity of the generated fractures. Thus, complex and organized geometries can be produced, and compared to existing mapped speleological networks as described in the cave and karst literature. Nevertheless, the presented analysis doesn’t show an appreciable influence of the fracture density on the karstogenesis evolution of the aquifer in terms of mean fracture aperture and total flow rate, at least for this setting. The role of the hydraulic conductivity of the carbonate matrix was also investigated. It can be noted that after a sufficient period of evolution (ranging from 1000 to 2500 years for a matrix hydraulic conductivity increasing from 10-6 to 10-4 m/s), the carbonate matrix exercises a decreasing control on the whole flow, while the enlarged fracture network becomes more and more influent. This fact appears clearly in the study of flow rate evolution. The

123

Carbonates Evaporites

more permeable the rocky matrix, the longer is its influence on the flow in the aquifer. Moreover, when the flow increases, this involves a faster enlargement of the fracture as the amount of young water supplied within the aquifer increases. These results show that ‘‘the porous continuum acts as a sink’’ in the case of a highly permeable matrix, which is consistent with former studies investigating the influence of different matrix hydraulic conductivities on the karstogenesis process (Bauer et al. 2000; Romanov et al. 2002; Kaufmann 2003; Kaufmann et al. 2010). Nevertheless, the choice of a time-stable matrix hydraulic conductivity prevents from simulating the evolution of the connectivity between fractures with time: indeed, the surrounding potentially fissured carbonate matrix remains soluble, and might evolve in parallel with the more developed fracture networks. This is a limitation for the realistic simulation of speleogenesis processes. Finally, this study shows that the model is able to simulate karstogenesis processes for a single stratum aquifer crossed by two orthogonal families of fractures and for laminar flow conditions. Further studies investigate simulations on multistrata setting, as the effect of various boundary conditions [namely diffuse and/or concentrated recharge conditions as presented by Lafare (2011)] and time varying boundary conditions. The question of the effect of different boundary conditions is indeed very important: in this study, the fixedhead gradient between input and output limit does not induce any limitation of the flow rate, which is not realistic. A different hydraulic boundary condition might therefore lead to different results, in terms of morphology and temporal evolution. A wider range of fracture networks should as well be performed, for example investigating the influence of different fracture orientations. Acknowledgements The authors thank the Conseil Re´gional of Languedoc Roussillon, as well as the BRGM and Carnot fund, for the financial support of this work carried out in the framework of the Antoine Lafare’s Ph.D. The authors would also like to thank Fabien Cornaton, presently Director, GMC at DHI-WASY GmbH (Germany), the CHYN of the University of Neuchaˆtel, (Switzerland), for permitting the using of the GroundWater modeling software, and for precious advices.

References Bakalowicz M (2006) Importance of regional study site conditions in elaborating concepts and approaches in karst science. Geol Soc Am Spec Pap 404:15–22. doi:10.1130/2006.2404(02) Bauer S (2003) Modeling of karst aquifer genesis: influence of exchange flow. Water Resour Res 39:1285–1297. doi:10.1029/ 2003WR002218 Bauer S, Liedl R, Sauter M (2000) Modelling of karst development considering conduit-matrix exchange flow. In: Proceedings of the ModelCARE’99 conference; calibration and reliability in groundwater modelling; coping with uncertainty, vol 265. IAHSAISH Publication, Zu¨rich, pp 10–15

123

Bloomfield JP, Barker JA (2005) MOPOD: a generic model of porosity development. Geol Soc Lond Spec Publ 249:73–77. doi:10.1144/GSL.SP.2005.249.01.06 Bloomfield JP, Barker JA, Robinson N (2005) Modeling fracture porosity development using simple growth laws. Ground Water 43:314–326 Cornaton FJ (2003) Deterministic models of groundwater age, life expectancy and transit time distributions in advective-dispersive systems (Ph.D. thesis). University of Neuchaˆtel, Neuchaˆtel, Switzerland Cornaton FJ (2007) Ground water (GW): a 3-D ground water flow, mass transport and heat transfer finite element simulator. University of Neuchaˆtel, Neuchaˆtel Cornaton F, Perrochet P (2006a) Groundwater age, life expectancy and transit time distributions in advective–dispersive systems; 2. Reservoir theory for sub-drainage basins. Adv Water Resour 29:1292–1305 Cornaton F, Perrochet P (2006b) Groundwater age, life expectancy and transit time distributions in advective–dispersive systems: 1. Generalized reservoir theory. Adv Water Resour 29:1267–1291 Cornaton FJ, Park Y-J, Normani SD, Sudicky EA, Sykes JF (2008) Use of groundwater lifetime expectancy for the performance assessment of a deep geologic waste repository: 1. Theory, illustrations, and implications. Water Resour Res 44:W04406. doi:10.1029/2007WR006208 Do¨rfliger N, Fleury P, Ladouche B (2009) Inverse modeling approach to allogenic karst system characterization. Ground Water 47:414–426. doi:10.1111/j.1745-6584.2008.00517.x Dreybrodt W (1996) Principles of early development of karst conduits under natural and man-made conditions revealed by mathematical analysis of numerical models. Water Resour Res 32:2923–2935 Dreybrodt W, Eisenlohr L (2000) Limestone dissolution rates in karst environments. In: Klimchouk A, Ford D, Palmer A, Dreybrodt W (eds) Speleogenesis: evolution of karst aquifers. National Speleological Society, Huntsville Dreybrodt W, Gabrovsek F, Romanov D (2005) Processes of speleogenesis: a modeling approach, Karst Research Institute, ZRC SAZU, Postojna, Slovenia Fleury P, Ladouche B, Conroux Y, Jourde H, Do¨rfliger N (2009) Modelling the hydrologic functions of a karst aquifer under active water management—the Lez spring. J Hydrol 365:235–243 Ford DC, Williams PW (2007) Karst hydrogeology and geomorphology. Wiley, Chichester Jeannin PY (1996) Structure et comportement hydraulique des aquife`res karstiques. Universite´ de Neuchaˆtel, Neuchaˆtel Josnin J-Y, Jourde H, Fenart P, Bidaux P (2002) A three-dimensional model to simulate joint networks in layered rocks. Can J Earth Sci 39:1443–1455 Jourde H (1999) Simulation d’essais de puits en milieu fracture´ a` partir d’un mode`le discret base´ sur des lois me´caniques de fracturation; validation sur sites experimentaux. Simulation of well testing in fractured media using a discrete model based on the mechanical laws of fractures; validation on experimental sites (Ph.D. thesis). Universite Montpellier II, Sciences et techniques du Languedoc, Montpellier, France Jourde H, Cornaton F, Pistre S, Bidaux P (2002) Flow behavior in a dual fracture network. J Hydrol 266:99–119 Jourde H, Fenart P, Vinches M, Pistre S, Vayssade B (2007) Relationship between the geometrical and structural properties of layered fractured rocks and their effective permeability tensor. A simulation study. J Hydrol 337:117–132 Kaufmann G (2003) A model comparison of karst aquifer evolution for different matrix-flow formulations. J Hydrol 283:281–289 Kaufmann G, Braun J (1999) Karst aquifer evolution in fractured rocks. Water Resour Res 35:3223–3238

Carbonates Evaporites Kaufmann G, Braun J (2000) Karst aquifer evolution in fractured, porous rocks. Water Resour Res 36:1381–1391 Kaufmann G, Romanov D, Hiller T (2010) Modeling three-dimensional karst aquifer evolution using different matrix-flow contributions. J Hydrol 388:241–250. doi:10.1016/j.jhydrol.2010.05. 001 Kiraly L, Jeannin P-Y, Sauter M (1998) Introducing to karst modelling; modelling in karst systems. Bull d’Hydrogeol 16:1–6 Klimchouk AB, Ford DC (2000) Types of karst and evolution of hydrogeologic setting. In: Klimchouk A, Ford D, Palmer A, Dreybrodt W (eds) Speleogenesis: evolution of karst aquifers. National Speleological Society, Huntsville Lafare A (2011) Mode´lisation mathe´matique de la spe´le´ogene`se: une approche hybride a` partir de re´seaux de fractures discrets et de simulations hydroge´ologiques (the`se de doctorat, Eaux Continentales et Socie´te´ - UM2). Universite´ de Montpellier 2, Montpellier, France Lafare A, Jourde H, Leonardi V, Pistre S, Do¨rfliger N (2009) Speleogenesis of Mediterranean karsts; a modelling approach based on realistic fracture networks. In: Hypogene speleogenesis and karst hydrogeology of artesian basins; proceedings, Ukrainian Institute of Speleology and Karstology Special Paper. Presented at the International conference; hyogene speleogenesis and karst hydrogeology of artesian basins, Ukrainian Institute of Speleology and Karstology, Simferopol, Ukraine (UKR), Chernivtsi, Ukraine, pp 75–81 Liedl R, Sauter M, Hueckinhaus D, Clemens T, Teutsch G (2003) Simulation of the development of karst aquifers using a coupled continuum pipe flow model. Water Resour Res 39:11

Mangin A (1975) Contribution a` l’e´tude hydrodynamique des aquife`res karstiques (Ph.D. thesis). Universite´ de Dijon, Dijon, France Mazzilli N, Guinot V, Jourde H (2012) Sensitivity analysis of conceptual model calibration to initialisation bias. Application to karst spring discharge models. Adv Water Resour 42:1–16. doi:10.1016/j.advwatres.2012.03.020 Palmer AN (1991) Origin and morphology of limestone caves. Geol Soc Am Bull 103:1–21 Palmer A (2000) Digital modeling of individual solution conduits. In: Klimchouk A, Ford D, Palmer A, Dreybrodt W (eds) Speleogenesis: evolution of karst aquifers. National Speleological Society, Huntsville Palmer AN (2001) Dynamics of cave development by allogenic water. Acta Carsologica 30(2):14–32 Romanov D, Dreybrodt W, Gabrovsˇek F (2002) Interaction of fracture and conduit flow in the evolution of karst aquifers. In: Martin JB, Wicks C, Sasowsky ID (eds) Proceedings of the symposium on karst aquifers: florida and related environments. Presented at the Symposium on Karst Aquifers: Florida and Related Environments, Karst Waters Institute, Charles Town, WV, pp 1–6 Sauter M (1992) Quantification and forecasting of regional groundwater flow and transport in a karst aquifer (Gallusquelle, Malm, SW. Germany) (Ph.D. thesis). Eberhard Karls Univ. Tu¨bingen, TGA, C13, Tu¨bingen, Germany Tritz S, Guinot V, Jourde H (2011) Modelling the behaviour of a karst system catchment using non-linear hysteretic conceptual model. J Hydrol 397:250–262. doi:10.1016/j.jhydrol.2010.12.001

123