openKARST: A novel open-source flow simulator for karst systems

TL;DR

The paper presents openKARST, a Python-based, open-source solver for fully transient conduit-flow in karst networks by solving the dynamic Saint-Venant equations for both free-surface and pressurized regimes. It combines a Preissmann slot mechanism with Churchill's friction formulation to enable smooth transitions across flow states and Reynolds regimes, while employing Crank-Nicholson and backward-Euler time integration with an adaptive timestep strategy. The approach is validated against analytical benchmarks and laboratory data, and demonstrated on a large-scale Ox Bel Ha cave network to illustrate scalability and the ability to capture rapid recharge and discharge dynamics. This tool advances karst hydrology by enabling high-resolution, transient simulations of complex conduit networks and lays groundwork for future coupling with matrix flow and solute transport analyses.

Abstract

We introduce the open-source Python-based code openKARST for flow in karst conduit networks. Flow and transport in complex karst systems remain a challenging area of hydrogeological research due to the heterogeneous nature of conduit networks. Flow regimes in these systems are highly dynamic, with transitions from free-surface to fully pressurized and laminar to turbulent flow conditions and Reynolds numbers often exceeding one million. These transitions can occur simultaneously within a network, depending on conduit roughness properties and diameter distributions. openKARST solves the transient dynamic wave equation using an iterative scheme and is optimized through an efficient vectorized structure. Transitions from free-surface to pressurized flows in smooth and rough circular conduits are realized via a Preissmann slot approach in combination with an implementation of the Darcy--Weisbach and Manning equations to compute friction losses. To mitigate numerical fluctuations commonly encountered in the Colebrook--White equation, the dynamic switching from laminar to turbulent flows is modeled with a continuous Churchill formulation for the friction factor computation. openKARST supports common boundary conditions encountered in karst systems, as and includes functionalities for network import, export and visualization. The code is verified via comparison against several analytical solutions and validated against a laboratory experiment. Finally, we demonstrate the application of openKARST by simulating a synthetic recharge event in one of the largest explored karst networks, the Ox Bel Ha system in Mexico.

Figures (11)

• Figure 1: Spatial discretization of the underlying network showing both network nodes (red circles) and computational nodes (black dots). Note that the spacing $\Delta x$ can be variable and different for each conduit of length $L$.
• Figure 2: Comparison of the numerical solution for the water depth at steady-state with the analytical solution Eq. \ref{['eq:analyical_delestre1']} for different resolutions $\Delta x=\SIlist[list-units=single]{1;5;10;25;50}{m}$. Notice that the channel base (lightgray) is only shown for a resolution of $\Delta x=1m$. Errors are shown for all resolutions with lighter gray for higher values of $\Delta x$. RMSE values are computed based on the percentage errors.
• Figure 3: Comparison of the numerical solution for the water depth at steady-state with the analytical solution Eq. \ref{['eq:analyical_delestre2']} for different resolutions $\Delta x=\SIlist[list-units=single]{1;5;10;50;100;200}{m}$. Notice that the channel base (lightgray) is only shown for a resolution of $\Delta x=1m$. Errors are shown for all resolutions with lighter gray for higher values of $\Delta x$. RMSE values are computed based on the percentage errors.
• Figure 4: Comparison of the numerical solution for different resolutions $\Delta x=$[list-units=single]1;5;10;25;50m for the water depth at steady-state with the analytical solution Eq. \ref{['eq:analyical_delestre1']}. A steady diffuse recharge of $r_0 = 0.001m s^{-1}$ is applied along the whole domain. Notice that the channel base (lightgray) is only shown for a resolution of $\Delta x=1m$. Errors are shown for all resolutions with lighter gray for higher values of $\Delta x$. RMSE values are computed based on the percentage errors.
• Figure 5: Comparison of steady-state discharge for three different roughness coefficients under pressurized flow conditions in a conduit of length $L=1000m$ with $\Delta x=100m$ (left) and corresponding errors (right). Circles represent the numerical solution obtained at an equilibrium at $t=2000s$ and lines the analytical solution computed with Eq. \ref{['eq:DW_corrected']}.