Karst conduit size distribution evolution using speleogenesis modelling

One of the critical aspects when modeling groundwater flow in karstic aquifers is to estimate the statistics of the size of the conduits, in conjunction with the connectivity of the karst conduit network. Statistical analysis can be performed on data gathered by speleologists, but a significant fraction of the karst conduit networks is not directly reachable, and therefore, the resulting statistics are incomplete. An alternative method to evaluate the inaccessible areas of a karst conduit network is to simulate numerically the speleogenesis processes. In this paper, we use a coupled reactive-transport model to simulate the evolution of a vertical section of a fractured carbonate aquifer and investigate how the statistical distribution of the fracture apertures evolves. The numerical results confirm that the karstification proceeds in different phases that were previously hypothesized and described (inception, gestation, development). These phases result in a multi-modal distribution of conduit aperture. Each mode has a roughly lognormal distribution and corresponds to a different phase of this evolution. These outcomes can help better characterize the statistical distribution of karst conduit apertures including the inaccessible part of the network.


Introduction
Karst aquifers develop when rock mass permeability is increased by dissolution processes (e.g., de Waele and Gutierrez 2022). The global behavior of these aquifers is strongly controlled by fast and mostly turbulent groundwater flow within karst conduits connected over large distances, and by slow recessions that can be fed by the fractured carbonate matrix. Modeling these processes is known to be difficult as highlighted, for example, by the Karst Modeling Challenge (Jeannin et al. 2021).
If we focus on physically based models, several numerical tools (such as MODFLOW-CFP, FEFLOW, or DISCO) can solve the flow and solute transport equations based on the knowledge of the karst conduit network geometry and the physical properties of the conduit and rock matrix (Reimann and Hill 2009;de Rooij et al. 2013;Kresic and Panday 2018). These tools require on one hand a 3D mesh representing the geometry of the karst conduit network as one-dimensional objects in a 3D space filled by matrix elements, and on the other hand the hydraulic properties of the matrix and karst conduits. While recent progresses have been made in developing methods allowing to generate plausible network geometries (e.g., Jaquet et al. 2004;Pardo-Igúzquiza et al. 2012;Rongier et al. 2014;Fandel et al. 2021;Banusch et al. 2022), the question of defining the proper statistical and spatial distribution of the conduit diameters remains much more open. These parameters are, however, critical to predict possible flow and solute transport in karst aquifers. The main source of data to constrain the conduit diameter statistics are direct measurements made by speleologists. The statistical and geostatistical analysis of these data (Pardo-Igúzquiza et al. 2012;Frantz et al. 2021) shows that the accessible conduits size has roughly a log-normal distribution and 360 Page 2 of 16 their size is spatially correlated along the conduit network. However, the data sets on which these studies are based are incomplete, they do not include conduit sizes that are too narrow to be explored by speleologists.
Speleogenesis modelling is a possibility to complement the data directly observable by speleologists and obtain some information about the statistical distribution of the smaller conduits. A closely related field is the study of the development of dissolution patterns in rough fractures or rock masses (Upadhyay et al. 2015;Lipar et al 2021). By using reactive transport models one can model the development of geological structures affecting the properties of fractures or the creation of karstic geomorphologies. Some particularly interesting results in this field were obtained by Upadhyay et al. (2015). They show that the evolving dissolution patterns are insensitive to the amplitude and correlation length of the initial aperture field. However, their study does not consider the transition from laminar to turbulent flow.
The speleogenesis modelling approach has been pioneered by authors like Dreybrodt (1990), Gabrovšek and Dreybrodt (2000) or Kaufmann (2003). It involves the coupled numerical simulation of groundwater flow and reactive transport that causes mineral dissolution in the host rock. The dissolution process alters the petrophysical properties of the rock, such as fracture apertures. There is a positive feedback mechanism, because the fractures enlarged by dissolution allow a larger flow rate. As the flow rate increases, the reactive water can penetrate deeper into the fracture network and accelerate the fracture growth rate. Dreybrodt et al. (2005) have simulated these coupled processes to gain insights into the evolution of conduits and synthetic fracture networks. Their simulations revealed the importance of 4 th order dissolution kinetics in the development of karst conduit. They explain the development of fractures hundreds or thousands of meters from the inlet of a fracture network. This type of approach has been applied at different scales and on different types of aquifers (e.g., confined, unconfined, coastal, deep, etc.) to better understand the effect of these various situations on the speleogenesis processes (e.g., Hubinger and Birk 2011;Perne et al. 2014;Kaufmann et al. 2014;de Rooij and Graham 2017;Cooper and Covington 2020). Among the known results, Hubinger and Birk (2011) show that the enlargement of fractures creates discrete pathways (connected network) with bimodal fracture aperture distributions, where only the largest fractures continue to grow after the breakthrough of a pathway connecting the inlet and outlet of the modelled network. The breakthrough is defined as the stage in fracture evolution when aperture has increased sufficiently to reduce friction loses and allow the flow regime to change from laminar to turbulent. The transition is followed by a significant increase in flow rate at the outlet. Hubinger and Birk (2011) also show that a unimodal fracture aperture distribution emerges if the recharge is severely reduced.
The aim of this paper is to investigate the statistical distribution and evolution of possible karst conduit sizes below what is accessible by a speleologist. For this purpose, we implement and use a speleogenesis model based on the concepts described in detail in Dreybrodt et al. (2005). The numerical model uses FEFLOW to solve the flow and transport equations. The code interacts with FEFLOW to couple the dissolution process with the reactive transport modeling. The initial stage is assumed to be a two-dimensional discrete fracture network to keep the computing time manageable. We test different initial fracture networks and recharge conditions to check the sensitivity of the results to the initial geometry and recharge. We then run the numerical model of speleogenesis to assess the evolution of the aperture of the fracture network, i.e., karst conduits. In particular we pursue the simulation after the breakthrough to investigate how the karst conduit aperture continues to evolve. As compared to previous studies, we show that the evolution of the system includes multiple phases that can be interpreted as phases of inception, gestation, and development as they were previously hypothesized by Filipponi (2009) but never yet described on the base of a numerical model. These multiple phases lead to multiple modes in the statistical distribution of the fracture apertures.
This paper is organized in three parts. We first present the conceptual model and numerical tool that was developed and used in this work. We then describe our numerical experiments and the results. Finally, we discuss the implications of those results and conclude.

Conceptual model and numerical approach
The numerical code employed for this study is an extension of a previous code developed by Maqueda (2017 The modelling approach implemented in KSP is based on the conceptual model presented by Dreybrodt et al. (2005). The speleogenesis model relies on three assumptions. The first assumption is that the walls of the fracture are assumed to be natural limestone (calcium carbonate) thus soluble by acidic water. The second assumption is the dissolution reaction occurs at the rock surface only, i.e., the effect of calcite dissolution is a retreat of the fracture wall. The third assumption is that the mineral dissolution is the only mechanism for fracture growth modelled. Other processes such as the erosion of the walls by suspended solids in water, or rock detachment along the conduit walls due to mechanical stress are not accounted for.
KSP uses the kinetics chemistry model of Dreybrodt et al. (2005). Calcite dissolution rate is expressed as a function of the ratio of dissolved calcium concentration C to calcite equilibrium concentration C eq estimated by applying equilibrium chemistry concepts (Appelo and Postma 2010). The calcium equilibrium concentration is a function of dissolved carbon dioxide (CO 2 ) concentration, temperature, and the presence of other dissolved mineral species. The lowest calcium concentration ratio (C/C eq ~ 0) yields the faster dissolution rate. A calcium concentration equal equilibrium concentration yields a dissolution rate equal to 0 mol/m 2 ⋅s. The main feature of the kinetics model is that the dissolution rate can be explained by a linear or 4th order reaction model depending on the C/C eq ratio. The model considers four dissolution rates for the following combinations of C/C eq ratio and flow regime. Equation 3 applies for laminar flow conditions when C/C eq < 0.9, while Eq. 4 applies when C/C eq > 0.9, Eq. 5 applies for turbulent flow conditions when C/C eq < 0.9, while Eq. 6 applies when C/C eq > 0.9: where R is the dissolution reaction rate [mol m −2 s −1 ], k l is the 1st order reaction rate constant, k n is the 4th order reaction rate constant, b is the fracture aperture [m], c is the calcium ion concentration [mol m −3 ], C eq is the calcium equilibrium concentration (mol m −3 ), and D is the calcium ion diffusion coefficient [m 2 s −1 ]. The mass transport is described by the advection-dispersion model applied to one dissolved species (calcium ion Ca 2+ ). The implementation of the mass transport model is presented in detail in FEFLOW documentation (Diersch 2013).
The change in fracture aperture caused by mineral dissolution was coded into KSP. KSP can simulate the transition between laminar and turbulent flow and can estimate fracture growth for both flow regimes. KSP switches from laminar to turbulent flow equation at a point in fracture evolution (cross section geometry, gradient), where both models yield the same resistance to flow. This is achieved by comparing in every fracture the flow velocity obtained by Hagen-Poiseuille model with the flow velocity computed by the Manning-Strickler model. When the turbulent solution yields a higher resistance to flow compared to the laminar solution for the same fracture geometry and gradient, the simulation changes the flow equation to turbulent for every specific fracture in the simulated domain. This approach offers a numerically stable solution for the transition.
The change in fracture geometry due to mineral dissolution (r = fracture wall retreat) is computed with Eq. (7), where R is the dissolution reaction rate depending on C/C eq ratio and flow regime, A is the fracture surface area, and V m is the molar volume of calcite rock: The numerical implementation of the method described above relies on the splitting of reactive transport and the wall retreat calculations based on a quasi-stationary state approximation (Litchner 1988). This approximation is the assumption that dissolution reaction rates are relatively slow; therefore, it is acceptable to extrapolate an estimated dissolution rate over a period of simulation time to compute fracture wall retreat. The approximation is implemented by simulating only reactive transport until solute concentrations in the fracture network tend to stabilize while using timesteps with a duration that ensures numerical stability. Only then, wall retreat is computed and KSP modifies the aperture of fractures in the FEFLOW model in a single time-step. Next, the reactive transport problem is left to run once again, until a new quasi-stationary state in mass concentration in the model is attained.

KSP benchmark test
To test the KSP code, we first compared the results obtained with KSP with those published by Dreybrodt et al. (2005). For this benchmark, the aim is to model the evolution of one single rectangular fracture having a width of 1 m, an initial aperture of 0.2 mm, and a length of 1 km. The hydraulic boundary conditions are: (i) a constant pressure head of 50 m at the inlet point, and (ii) a constant head of 0 m at the outlet. The solute boundary condition is water fully unsaturated with calcite at the inlet. The calcite equilibrium concentration is 2 mmol/liter. The linear and 4 th order reaction kinetics constant values are 4⋅10 -7 and 4⋅10 -4, respectively. Before comparing the results, note that there are two fundamental differences between the Dreybrodt et al. (2005) simulations and our implementation. In Dreybrodt et al. (2005), the flow occurs only in the fracture and the hydrodynamic dispersion of the solute is not accounted for. Whereas in our simulation, the fracture is embedded in a porous rock matrix with low hydraulic conductivity (< 10 -6 m/s), where flow and solute transport still occur at a minimal rate, and hydrodynamic dispersion, representing the heterogeneity in water flow velocities in the fractures, is accounted for within the advective-dispersion formulation in FEFLOW. Figure 1a presents the results of the evolution of fracture aperture in both the KSP simulation (dashed lines) and the Dreybrodt et al. (2005) simulation (solid lines). At simulation times of 13,100 and 17,800 years, the fracture aperture in both simulations is nearly identical. The transition to turbulent flow occurs in both simulation at around 18,850 years, whereas the fracture apertures in the KSP simulation are slightly smaller than the aperture at the benchmark model (blue lines). The largest difference in fracture aperture is observed at a simulation time of 19,032 years. We interpret the deviations as resulting from the differences in implementation due to the presence of porous media in KSP and the different transport equations. At a simulation time of 19,152 yrs. the difference in fracture aperture is reduced after flow becomes turbulent. Figure 1b presents the evolution of flow rate for both KSP simulation and Dreybrodt et al. (2005) simulation. The simulated flow rates with KSP are comparable to the benchmark before and after the transition to laminar flow (nearly vertical increase in flow rate). Only a small difference is observed by the end of the simulation probably, because flow and solute transport in porous media in our simulation. The conclusion from the benchmark test of the KSP is that the code reproduces the main trends of fracture aperture and is able to model the transition between laminar and turbulent flow to reproduce the flow rate increase after breakthrough Dreybrodt et al. (2005) simulation. Finally, we consider that this test shows that KSP is capable to simulate reasonably well-fracture aperture evolution from 10 -4 m to 10 m.

Model setup
KSP was used to investigate the evolution of a simplified karst aquifer under phreatic conditions. The geometry is a rectangular vertical cross section through a fractured carbonate formation (Fig. 2). Recharge is acting on the top of the system, and a spring is located on the lower part of the right side of the domain.
The 2D model domain has a size of 2 km horizontally and 500 m vertically. Initial fracture networks 1 and 2 within the model domain were generated using a simple object-based model that we implemented in python. Two different fracture networks were considered (Table 1). For each one, 4 families of discontinuities (fracture sets) were defined. The sub-horizontal family represents bedding planes or sub-horizontal tectonic discontinuities. The sub-vertical family represents sub-vertical tectonic discontinuities. In addition, there are two families of conjugated fractures. All the fractures are simulated independently. Their position is generated following a Poisson random point process with a density that is different for every fracture family. The distribution of the length of the fractures follows a truncated power-law distribution, with an exponent that has been kept constant. The orientations follow a von Mises distribution for each fracture family. All the parameters of those statistical distributions are provided in Table 1.
In total, the fracture network 1 has 6,257 fractures with a total accumulated fracture length of 57,628 m; the fracture network 2 has 5,554 fractures with a total accumulated fracture length of 39,727 m. The differences between fracture networks 1 and 2 are the density and minimum length of fractures (Table 1). This difference is well-visible on Fig. 3,  where fracture network 1 looks denser and with shorter fractures when compared to fracture network 2. Fracture apertures have been shown to be variable in space and spatially correlated at multiple scales (Tatone and Grasselli 2012). To generate a simple but plausible initial distribution of apertures, we used the sequential gaussian simulation method within the SGEMS software (Remy et al. 2009) to generate the logarithm of the apertures as a random multi-Gaussian field, with a Gaussian variogram model with a range equal to 150 m in the horizontal direction and 50 m in the vertical direction. The random multi-Gaussian field is created as a raster file which is used to assign initial apertures to the fracture networks in the numerical model The final distribution of the aperture has a log-normal distribution with a mean of 1.68⋅10 -4 m and a variance of 3.8⋅10 -9 m 2 ; the initial apertures have values ranging between 5.0•10 -5 m and 5.0•10 -4 m (Fig. 3). The fractures are discretized into smaller elements, whereas each fracture element has a unique initial aperture value to improve the convergence and numerical accuracy of the numerical solution.
When the flow becomes turbulent, we also need to define the value of the Strickler coefficient. For this purpose, we use a range of values previously identified in the scientific literature. Jeannin (2001) estimated a Strickler coefficient of 20 based on karst conduit geometry and flow measurements for the Hölloch cave system in Switzerland. An estimation of the Strickler coefficient was also done for the Devil's Icebox-Connor's Cave system in the USA based on flow and geometry data and yields a Strickler coefficient between 28 and 150 (Peterson and Wicks 2006). Perne et al. (2014) used values between 50 and 100 for simulations of the transition from pressurized flow to free-surface flow in karst conduit networks. Based on the aforementioned references, we assume a Strickler coefficient of 50 which should be representative of karst conduit networks.
For the flow boundary conditions (Fig. 2), we assume that the amount of precipitation is constant. In the initial state of the system, the amount of precipitation is higher than what can enter in the fractures, because the initial fracture permeabilities are low. Therefore, the total amount of available recharge by precipitation cannot completely enter the system and a large proportion of it is eliminated by surface runoff. This process maintains the fractures fully saturated.
For each fracture network we considered three simulation scenarios with different values for the constant pressure head boundary condition: case (a) h = 33 m, case (b) h = 100 m and case (c) h = 300 m. The different values of constant pressure head will allow to assess its effect on the evolution of the fracture networks. When the system evolves, the fractures are enlarged by dissolution, the network becomes much more permeable, but the recharge of the system cannot become larger than the rate of precipitation on the upper surface. Therefore, we implemented a constraint of maximum inflow rate in addition to the constant head boundaries. The maximum flow rate is set to 1 l/s per recharge node. This hydraulic boundary condition is implemented on the nodes,  Fig. 3.
The outlet of the system (spring) is located at the intersection between one fracture and the right boundary near the bottom of the model domain (Fig. 3). The spring is modeled using a prescribed constant head h = 0 m. The location of the spring is also represented as a blue circle on the right side on the bottom of the domain in Fig. 3.
Mineral dissolution (reactive transport) is simulated with a variable time-step duration in the order of hours controlled by FEFLOW solver to ensure a stable numerical solution. The solute concentration arrives to a quasi-stationary condition (stable concentration) after 50-100 days of simulation. The simulation time-step duration for the simulation of fracture enlargement is fixed at 7,300 days (~ 20 years). The total simulated period of fracture growth is between 5 and 10 K years depending on simulated domain size, fracture geometry and hydraulic boundary conditions. Figure 4 presents the fracture aperture distributions for all the simulated cases at the simulation time when the maximum flow rate of 1 l/s per recharge node is attained. Even if the initial fracture aperture is identical for fracture networks 1 or 2, the final structure of the enlarged fracture network depends on the hydraulic conditions and can favorize different fractures to grow faster. With smaller hydraulic gradient (e.g., case 1a), the resulting enlarged fracture networks tend to be "simpler", i.e., it contains less conduits (enlarged fractures) compared to resulting networks with stronger hydraulic gradients (e.g., cases 1b, 1c). The same outcome is observed for scenarios 2a, 2b, and 2c, where stronger gradient causes a more "complex" network of turbulent flow fractures. Figure 4 also shows that karst aquifers do not need to develop all over a karstified rock mass but only in limited areas. For example, larger fractures or karst conduits did not develop significantly below the outlet elevation. Karst conduits are also absent in areas distant from the outlet. In these areas without karst conduits, the hydraulic behavior is determined by flow along fractures (fracture aquifers), while in those with karst conduits, these dominate the aquifer hydraulics (karst aquifers).

Fracture aperture evolution
The networks illustrated in Fig. 4 are the results of multiple discrete breakthroughs corresponding to the connection of several new clusters or parts of the fracture network to the turbulent and rapid flow path system between the recharge area and the outlet. This evolution is illustrated in Fig. 5.   Fig. 4 Comparison of the final geometry (developed length) of the fracture network for all simulation cases (a-c) for fracture networks 1 and 2 The figure shows the evolution of fracture apertures for case 1a at four important and distinct simulation times. At time T = 1560 years, a small cluster of fractures located above the outlet has just and rapidly transitioned from laminar to turbulent flow, their apertures became much larger than the initial values (red colors). The histograms only present aperture distribution, the relationship to the transition from laminar to turbulent flow is presented in the videos included as supplementary material for this article. These fractures connect the outlet to the recharge area. At T = 1855 years, a new cluster of fractures has just made the same transition from laminar to turbulent flow (fractures in red within the highlighted dashed rectangle) and did connect to the first cluster. Later, at T = 3618 years and T = 5667 years, the same process repeated itself with the connection of additional clusters of fractures, highlighted again by the dashed rectangles in Fig. 5. This process involving the sudden connection of novel parts of a connected network of highly open fractures can be fully observed on the video animations provided as electronic supplementary information. These animations also permit to see that the same behavior occurs (with some variations) in all the simulated cases and include the isolines of hydraulic head.
The breakthrough of enlarged fractures influences the fracture aperture distribution as shown on the histograms of fracture aperture (Fig. 6). Frequency is shown in the vertical log scale, while the fracture aperture in meters is presented in horizontal log scale. For the case 1a, the initial fracture apertures (T = 0 years) have a log-normal distribution with a mode or peak indicated with the 0 symbol on the top of the figure. At T = 1,560 years, the 1st breakthrough causes an aperture mode (peak 1) of ~ 0.1 m. At T = 1,855 years, the 2nd breakthrough causes a 2nd aperture mode (peak 2) of ~ 0.1 m. By then, the fractures of the 1st breakthrough (peak 1) have grown and have a new mode of ~ 0.3 m. At T = 3,618 years, the third breakthrough causes a new aperture mode (peak 3) of ~ 0.1 m. At this time, peaks 1 and 2 have merged into a single peak of mode ~ 1 m. At T = 5,667 years, the 4th and final breakthrough occurs and yet another mode (peak 4) of ~ 0.2 m emerges. By then, peaks 1, 2, and 3 have nearly converged to an aperture between 1 and 3 m and can be regarded as karst conduits. At T = 5,667 years, the inflow rate stabilizes (see Fig. 4 case 1a), and the simulation is run until T = 7,115 years. At this evolution stage, the 4 th breakthrough fractures (peak 4) almost converge with previous peaks 1, 2, and 3, and a 5th mode emerges with an aperture ~ 0.2 m emerges. Since the inflow rate has stabilized, the flow in the background fractures does not transition from laminar to turbulent flow and it is not expected that fractures of peak 0 grow into conduits with turbulent flow, i.e., into karst conduits. Figure 7 presents the distribution of the final fracture aperture for all the simulated cases. Log scales are the same as Fig. 6. A multimodal distribution (several peaks) is observed in the fracture aperture distributions. The various modes are the result of the same evolution process presented for case 1a in Fig. 6. The multimodal distribution is the result of every new phase of karstification when a new group of enlarged fractures are connecting to the cluster of fractures previously connected to the spring. When comparing the outcomes of the six simulation cases a general trend of three modes of fracture aperture emerges: • A peak composed by fractures with apertures larger than 1 m (red color in Fig. 7). This peak represents all the fractures where flow transitioned from laminar to turbulent flow. These fractures keep growing after the stabilization of flow rate which causes this peak to be detached from the rest of the distribution. This peak represents fractures that developed to caves, i.e., these are the ones that are large enough to be explored in nature by speleologists. • A second peak is observed for fractures in the range of around 0.1-0.3 m (orange color in Fig. 7). These fractures are connected to the larger fractures of the first peak. However, they stopped growing once the flow rate stabilized, because the recharge of reactive water is drained by larger fractures (peak 1) to the spring. • Finally, a third aperture mode still in the range of the initial fracture aperture between 5.0•10 -5 to 5.0•10 -4 m (blue color in Fig. 7). The flow regime in these sub-millimeter fractures remains laminar. Figure 8 presents the evolution of the outflow rate at the spring for all the simulated cases. This evolution can be divided into three main phases:

Flow rate evolution
• An initial phase with a slow exponential increase before any connected fracture network has been sufficiently enlarged to transition from laminar to turbulent flow. During this initial phase, the flow rate can be approximated by the product of the total head gradient and an effective hydraulic conductivity that depends on the initial structure of the fracture network and the distribution of the apertures. The differences in flow rate between the cases are linearly proportional to the head gradient. The second phase occurs when fracture flow transitions from laminar to turbulent flow. This is characterized by a sudden increase in flow rate. At this stage, a complete flow path became turbulent and connects a part of the recharge area with the outlet. During that phase, we observe multiple successive breakthroughs for scenarios 1a, 1b, 1c and 2a corresponding to the sequential connection of different areas of the fracture networks with the outlet, which causes a stepwise increase of the flow rate at the outlet (see Fig. 5). These flow rate steps are less visible in cases 2b and 2c (Fig. 8). • During the third phase, the flow rate stabilizes and reaches the maximum amount of recharge available. Further enlargement of the fracture aperture is minor and has little impact on the flow paths. At this stage the hydraulic gradient within the aquifer area with the turbulent flow path tends progressively toward zero. Table 2 presents a summary of the values of the flow rates at the outlet fracture at the 1st breakthrough time and the time of stabilized flow rate. As expected, stronger pressure heads yield in a shorter breakthrough time (case 1c < 1b < 1a). Network 1, in spite of being more complex than network 2, displays a shorter time for flow stabilization, because the maximum recharge flow rate is greater than for network 1 (46 > 26 l/s).

Network statistics
To compare the complexity of the resulting turbulent flow fracture networks, we computed the following summary statistics (Table 3): • number of fractures experiencing turbulent flow and ratio of this number to the total number of fractures in the network. • sum of the length of fractures experiencing turbulent flow and ratio of this length to the total fracture network length. Fig. 8 Flow rate evolution and breakthrough times for all simulation cases. The multiple breakthroughs presented in Fig. 5a for case 1a are highlighted with arrows. The first increase in flow rate is the consequence of the 1st breakthrough (Fig. 5a), 2nd flow rate increase is caused by second breakthrough (Fig. 5b), and so on until the last increase in outflow rate caused by the 4th and final breakthrough as presented in Fig. 5b These summary statistics were estimated at two simulation times: (i) at the stabilization of the total outflow flow rate (in Table 3 presented in bold font), and (ii) at the simulation end (in Table 3 presented in normal font). For all the simulated cases, only between 13% and 20% of all the fractures transitioned from laminar to turbulent flow in the simulated fracture network and only a fraction of those grew into karst conduits. Whereas at both fracture networks 1 and 2, the ratio of fractures in turbulent flow is higher for the cases with higher hydraulic gradient. It is noted that in cases 1a, 2a, 1b and 1c, the number of fractures, where flow becomes turbulent after the stabilization of flow rate is between 0.25% and 1.7%. Whereas in cases 2b and 2c the number of fractures in turbulent flow is the same at the flow stabilization time and final simulation time. This can be understood as follows: when the maximum flow rate is achieved, the pressure head required to push reactive water through the smallest fractures in the network decreases and approaches zero. Then, the system does not evolve anymore significantly.
The second statistic in Table 3, the ratio of the sum of the length of the turbulent fractures shows that the turbulent fractures represent between 13% and 20% of total fracture network length. In this aspect, the cases with higher pressure head are the cases with more turbulent fractures. Notably, the cumulated length of fractures in turbulent flow is the nearly the same at flow stabilization and final simulation time demonstrating that once the pressure head approaches to zero, the development of new paths between inlets and outlets stops.
The relation of the evolution of the turbulent fracture network and the hydraulic gradient is not straightforward and does not follow simple intuition. For example, some flow paths that were enlarged with low initial hydraulic gradient have not been enlarged with higher hydraulic gradients. Depending on the hydraulic gradient, different fractures were enlarged and reached turbulent flow conditions (Fig. 4). This phenomenon is likely the result of the non-linearity of the underlying physics and chemistry. If the flow would only follow a linear equation, such as Poiseuille, increasing the initial hydraulic gradient will increase the flux in all fractures in a proportional manner. Because of the transition from laminar to turbulent flow, the resistivity, and the behavior of the enlarged fractures, changes in a non-linear manner. This effect is then amplified by the dissolution process and by the changes in boundary conditions on the top of the domain when the flow rate becomes higher than the maximum recharge. All these interactions between the processes are explaining the complexity of the observed behavior.

Conclusions and discussion
In this paper, we simulated the evolution of the aperture of a fracture system using a new implementation of the speleogenesis model of Dreybrodt et al. (2005). This model was implemented as a plugin within FEFLOW. As compared to the original model, the new implementation differs on some aspects related to the solute transport equation that includes a dispersion term in FEFLOW. There are also some differences in the way that the transition between laminar and turbulent flow are implemented. It would be useful to run several additional experiments to analyze in detail the impact of the porous matrix and dispersion equation on the results in the future to clarify the mechanisms controlling those differences. Nevertheless, the two models give rather similar results in the case of a single fracture. We then use this model to study how two 2D discrete networks of fractures evolve over time under different hydraulic boundary conditions.
One of the major results of this study is that we observe that the karstification process does not occur homogeneously and regularly in the catchment, but instead it proceeds in a series of multiple breakthroughs corresponding to different phases of karstification. These results illustrate for the first time in a quantitative manner a process that was proposed earlier by Lowe (1992) and extended by Filipponi (2009). The original concept included 3 phases in cave development: (i) inception, (ii) gestation, and (iii) development, as illustrated in Fig. 9. The inception phase was described as the start of dissolution in the fractures under phreatic conditions. The outlet or spring was assumed to be the consequence of a valley incision in a soluble rock massif. The outlet was hypothesized to organize the groundwater flow and at time 1, a cave gestation zone was assumed to emerge. This is very similar to what we observe during the initial phase of our simulations. At time 2, the breakthrough occurs, and the new karst conduit network acts as a "spring area" for the upstream section of the model which is comparable to the 1st breakthrough in our simulation. The karst conduits (cave development phase) offer less resistance to flow, thus the water table drops (observed as pressure head decrease in our simulation), and the gestation and inception zones move upstream. At time 3, the cave development keeps advancing upstream, which is comparable to the 2nd, 3rd, and 4th breakthroughs in our simulation cases 1a, 1b, 1c and 2a. In our simulations, the variability of the distances between the inlets and the outlets of the fracture networks resulted in several stages of breakthrough time when different regions of the network did connect the inlets and outlets. In our simulations the fractures are always full saturated.
The second important result that we obtained is that the distribution of fracture apertures is generally multimodal. Depending on the flow conditions and initial geometry, the distribution can be dominated by two main modes or more. This result is comparable and generalizes the results obtained by Hubinger and Birk (2011) who conducted numerical simulations of the growth of fracture networks by mineral dissolution. They found that during the initial phase of their simulations there was only 1 mode of fracture aperture. At later simulation times the modes representing enlarged fractures merge. In our simulations, the initial lognormal distribution of fracture size with one mode evolves into a log-normal mainly tri-modal distribution. The first mode can include several sub-families, but it generally corresponds to the fractures with an aperture larger than 1 m with flow in turbulent regime and which can be regarded as karst conduits. The second mode represents fractures that intersect the main conduits but where flow after the stabilization of flow rate is not sufficient to drive reactive water into them, and therefore, the enlargement of these fractures stopped. The third mode represents the initial fractures which experienced minimal growth and stay in laminar flow conditions during the entire simulated period.
As we discussed in the previous paragraph, the different modes have an essentially log-normal distribution of fracture aperture. This is similar to what has been observed when computing statistics of explored conduits (Maqueda 2017;Frantz et al. 2021). We also observe in our results that there is no trend of larger conduits found upstream or downstream in the network. This is again consistent with the analysis of field data by Frantz et al. (2021). However, the spatial distribution of the conduits radius in natural caves is known to be spatially correlated along the conduit network (Pardo-Igúzquiza et al. 2012;Frantz et al. 2021). Some branches of the conduit networks have larger conduit sizes, while others are narrower. A geostatistical analysis of the fracture Fig. 9 Conceptual model for the development of karst conduit network from Filipponi (2009) apertures obtained in our simulation has not been conducted, by a visual inspection of the results (e.g., in Fig. 5) show that there is little variability of the apertures along the cluster of connected enlarged fractures. Our model does not seem to be capable of reproducing the variability of these apertures even if we started from a random aperture distribution. This is partly a consequence of the chosen model of mineral dissolution kinetics during laminar flow conditions (before breakthrough time). Under laminar flow conditions, the dissolution rate is controlled by mass diffusion as seen in Eq. 3 (Dreybrodt et al 2005). The larger the fracture, the slower the dissolution reaction. Therefore, up to the breakthrough time and transition to turbulent flow, smaller fractures grow faster than larger fractures. However, after the stabilization of the flow rate, the reactive water flows mainly in the karst conduits (large fracture) and the relatively faster dissolution reaction kinetics becomes irrelevant, because no reactive water flows through them. Therefore, they stop growing. Other reasons possibly explaining that the model produces conduits of relatively homogeneous diameters are: (i) that the rock is assumed to be completely homogeneous, and (ii) the model only considers mineral dissolution, while in reality, mechanical effects also contribute to fracture growth.
Among the other results obtained in this study, we observe that the breakthrough occurs at a flow rate of similar order of magnitude for all the simulated cases. This behavior is interpreted as a result of the flow rate, dissolution reaction kinetics and fracture initial aperture interacting with each other. When fractures grow and flow velocity increases, reactive water can penetrate deeper into the fracture network. This penetration distance is what drives the breakthrough. Since breakthrough occurs at approximately the same flow rate, regardless of initial hydraulic boundary conditions, we conclude that breakthrough flow rate is a characteristic of the initial fracture network, mainly the resistance to flow of the fracture network. And finally, significant differences were observed for the simulations with same initial fracture network geometry and aperture but different initial hydraulic gradient. We quantified the complexity of the resulting karst conduit network in terms of number of fractures experiencing turbulent flow and concluded that when a larger hydraulic gradient push more water into more small fractures before the 1 st breakthrough time, a more complex conduit network tends to emerge.
Before concluding, we are aware that the results and conclusions expressed in this paper were obtained only for a few 2D models. To get more confidence in our results, it would be useful to run an ensemble of stochastic simulations and repeat a similar analysis on a large number of models. Furthermore, since the dimensionality of a model affects strongly the connectivity of a network, we expect that some of our conclusions may not remain valid in 3D. Exploring these effects would be rather straightforward but would imply much larger computation time and this was not possible in the framework of the present study.
Finally, the results presented in this paper should help understand how the statistical distribution of fracture apertures or conduit diameters can evolve during the speleogenesis. These results can serve to guide the selection of karst conduit aperture distributions that one could use for modeling groundwater flow and solute transport in karstic aquifers using physically based models. These results can have practical implications for diverse applications, including, for example, environmental impact assessments or inflow risk analysis when planning the construction of underground infrastructures (e.g., Casagrande et al. 2005;Filipponi et al. 2012).