Through a research initiative between Waterloo Hydrogeologic and University of Neuchatel, Switzerland, a new methodology to model karst conduits was developed in a stochastic framework. The research work has concluded, and now we would like to hear from you on the relevance of applying this methodology to groundwater modeling codes.

Introduction
Karst aquifers are important groundwater resources and are considered as the main freshwater resource for 25% of the world’s population. Karst aquifers are very heterogeneous. The dissolution of calcite by infiltrating water produces preferential flow paths (conduit networks) within these aquifers. The presence of these large voids can present geotechnical risks (tunnels, dams, etc) and create flow paths which make karstic springs particularly vulnerable to contamination. Quantifying the uncertainty is therefore essential, but the available modeling techniques for karst aquifers make it challenging to allow this kind of studies. Stochastic approaches have proven their efficiency in uncertainty estimation in many other fields of application. The present thesis proposes a new methodology to model karst aquifers, in a stochastic framework. The resulting models are consistent with the actual knowledge about these aquifers and their particular genesis processes (speleogenesis).
Methodology
The proposed methodology is divided into several steps:
- First, a geological model is built to model the large scale geometry of the aquifer.
- Secondly, the initial heterogeneities (inception horizons) of the aquifer formation are modeled stochastically, because they influence the small scale geometry of the conduits network.
- The stochastic conduit networks are then simulated using the physical assumption that water follows a minimum effort path while minimizing the total energy of a system.
- Finally, these models are used to simulate groundwater flow and transport. The models can be subsequently used in an inverse framework to estimate the uncertainty.

1. Geological Model
The construction of a realistic geological model is used to provide the information about the 3D localization of karstifiable formations and their orientation.

Figure 1a) and 1b): Geological model of the Valley of les Ponts-de-Martel
2. Discrete Fracture Network (DFN)
The geological model is used to generate stochastic Discrete Fracture Networks (DFN) with a new approach: the theoretical families of fractures that are expected in a folded environment are modeled with the consistency of the local geological orientation. The fractures (simple rectangular planar objects) are first modeled in an unfolded environment, and then placed in the geological model and rotated according to the local orientation of the geology. The fracture number and size are distributed according to a power law as defined in the Universal Fracture Model (UFM) for the dense regime.

3. Stochastic Conduits Network Modeling
The stochastic conduit networks are modeled using a pseudo genetic methodology. After the geological model is defined, and the initial heterogeneities are built, a study of the regional hydrology / hydrogeology is conducted to identify the potential inlets and outlets of the system, the base levels and the possibility of having different phases of karstification. The conduits are then generated in an iterative manner using a Fast Marching Algorithm. A probabilistic model can be used to represent the degree of knowledge available and the remaining uncertainty depending on the available data. The conduits are assumed to follow minimum effort paths in the heterogeneous medium from sinkholes or dolines toward springs. The search of the shortest path is performed using a Fast Marching Algorithm (FMA). This process can be iterative, allowing to account for the presence of already simulated conduits and to produce a hierarchical network. The final result is a stochastic ensemble of 3D karst reservoir models that are all constrained by the regional geology, the local heterogeneities, and the regional hydrogeological/hydrological conditions.
Several levels of uncertainty can be considered (large scale geological structures, local heterogeneity position of possible inlets and outlets, phases of karstification). Compared to other techniques, this method is fast, accounts for the main factors controlling the 3D geometry of the network, and can be conditioned to available field observations.

of FMA from the spring location (green point) to every point of the domain.
The black lines are the first conduits generated following the time map;
color: blue (low values) to red (high values). c) Third iteration of the algorithm.
Figure 2a) Initial fracture network, b) first iteration of the algorithm, initial time map resulting
of FMA from the spring location (green point) to every point of the domain.
The black lines are the first conduits generated following the time map;
color: blue (low values) to red (high values). c) Third iteration of the algorithm.


the vadose zone are generated toward a base level first, (b) and then the saturated conduits
are simulated toward the spring (b). The dotted line represents the phreatic surface.
The filled circle represent the spring, the grey lines represent the fractures and the bedding planes.
Figure 3a) result of one 2D cross-section simulation: (a) the conduits (black lines) within
the vadose zone are generated toward a base level first, (b) and then the saturated conduits
are simulated toward the spring (b). The dotted line represents the phreatic surface.
The filled circle represent the spring, the grey lines represent the fractures and the bedding planes.
4. Groundwater Flow and Transport Simulations
The karst conduits networks are then used to simulate groundwater flow and transport.
While the initial work was carried out using the finite element code GROUNDWATER, allowing for turbulent flow simulations in the conduits, recent applications have used the MODFLOW simulator. The realistic karst models (geometry and parameters) are converted to equivalent porous media and translated to MODFLOW packages (grid, properties, and basic package). The figures below show the corresponding MODFLOW properties for the conduits, as generated by the stochastic karst generator.

2D View of MODFLOW properties representing the conduit networks generated using SKS

3D View of MODFLOW properties representing the conduit networks generated using SKS
A lumped parameter model is used to model the recharge function. Specified head boundary conditions are used for the springs, specified flux boundary conditions for the point inlets of the aquifer (dolines, sinkholes), and specified flux for the diffuse recharge. The diffuse recharge represents only a small percentage of the total recharge. The point recharge is distributed over all the inlets depending on their catchment area. Transport transient simulations are used to simulate the tracer tests. A standard workflow is used: first the heads from a steady-state flow simulation are used as initial heads for the transient simulation. Secondly, several transport simulations are run (their number depend on the number of tracer tests that have to be simulated). The duration is compared to the transient flow simulation (only several weeks, or months, depending on the rear tracer tests results).
The flow and transport simulations and the geometrical realizations can be combined in an inverse framework. A possibility is to produce many (hundreds) geometrical realizations (conduits networks), and then to optimize the flow parameters of each one and study the uncertainty on this base. In this work, numerical tests have been done to estimate the feasibility of this framework. The results of this test show that including the transport results in the computation of the objective function can consistently improve the inverse results. Tests have also been run using PEST (Parameter Estimation) to retrieve realistic physical parameters of a realization. The results shows that this algorithm can be used in this case, and that it allows to save a considerable amount of run time in comparison with systematic search.
For more information
Please contact us if you are interested in discussing how this methodology can be extended to MODFLOW packages, such as MODFLOW-CFP (Conduit Flow Package) or CLN (Connected Linear Network package, available in MODFLOW-USG). Also send us any general questions or comments: sales@waterloohydrogeologic.com.
References
Andrea Borghi, Philippe Renard, Sandra Jenni. A pseudo-genetic stochastic model to generate karstic networks. Journal of Hydrology 414–415 (2012) 516–529.
C. Vuilleumier, A. Borghi, P. Renard, D. Ottowitz, A. Schiller, R. Supper, F. Cornaton. A method for the stochastic modeling of karstic systems accounting for geophysical data: an example of application in the region of Tulum, Yucatan Peninsula (Mexico). Hydrogeology Journal, 2012. DOI 10.1007/s10040-012-0944-1.
How to model realistic 3D karst reservoirs using a pseudo-genetic methodology – example of two case studies. A. Borghi, P. Renard, S. Jenni. ISKA 2010: 4th International Symposium on Karst – Málaga, 27th-30th April 2010
A pseudo-genetic stochastic methodology to model realistic karst reservoirs. Andrea Borghi, Philippe Renard, Sandra Jenni. IAHR 2010: International Groundwater Symposium – Valencia, September 2010
Fabien Cornaton, 2007, “GROUNDWATER: A 3-D Ground Water and Surface Water Flow, Mass Transport and Heat Transfer Finite Element Simulator”, University of Neuchâtel