Simulating precipitation-induced karst-stream interactions using a coupled Darcy–Brinkman–Stokes model

The interaction mechanism between karst aquifers and streams remains unclear, particularly regarding the impact of dynamic groundwater saturation processes under variable precipitation. This challenge hinders the accurate modeling of karst hydrology. This study developed a Darcy–Brinkman–Stokes model to analyze these complex interactions. The model integrates water-air two-phase flow and employs multiple water retention models to characterize variably saturated flow in porous media. We validated the DBS approach by comparing its numerical results against the MODFLOW-Conduit Flow Process v2 for generalized karst models. The key conclusions are as follows:

• Rainfall intensity is the dominant driver of the interaction. Higher intensities lead to more complex processes, involving multi-media collaborative recharge and shifting discharge contribution ratios from different media.

• During consecutive rainfall events, groundwater stored in porous media (matrix) significantly influences subsequent stream levels, whereas conduit storage shows negligible carry-over impact due to rapid drainage.

• Uncertainty analysis demonstrated that conduit geometry, epikarst permeability, and matrix porosity differentially influence system hydrology, controlling the magnitude, timing, and distribution of peak discharges.

The validated DBS model is a robust tool that accurately depicts the complex two-phase interactive flows (including infiltration, overflow, and recession) controlled by dynamic saturation. It successfully reveals the dynamic interactions between the epikarst, conduits, matrix, and stream, which is essential for understanding and managing karst water resources.

1 Introduction

Karst aquifer is not only a repository of substantial freshwater resources (Li et al., 2017; Ford and Williams, 2007; Sivelle et al., 2021), but also provides drinking water for 10 % to 25 % of the global population (Longenecker et al., 2017; Goldscheider et al., 2020; Mahler et al., 2021). However, karst-developed areas feature intricate pore structures and fractures (Kuniansky, 2016), leading to pronounced heterogeneity and anisotropy in the movement and storage of water within them (Zhang et al., 2020). In particular, the complex coupled flow involving various flow paths such as karst conduits, sinkholes, and epikarst, along with porous media, further intensifies the nonlinear recharge and discharge processes and the formation of preferential flow paths in the karst aquifer. With seasonal variations in precipitation intensity, the heterogeneity of the groundwater flow field is further exacerbated, and water levels in the karst aquifer and stream fluctuate, leading to complex interactions between the aquifer and stream (Bonacci, 2015). Unveiling the interaction mechanism between the karst aquifer and stream under varying precipitation intensities is crucial for assessing the storage of water resources in karst regions (Gao et al., 2021; Guo and Jiang, 2020).

The interaction process between the karst aquifer and stream is significantly influenced by karst media. In epikarst where the soil layer is shallow and dissolution weathering is pronounced, most precipitation can directly recharge the karst aquifer (Lee and Krothe, 2001; Okello et al., 2018). Karst conduits and sinkholes are important media involved in karst hydrological cycle. Together, they form a complex network for groundwater recharge and drainage. Surface water collected into sinkholes can directly recharge the karst aquifer (Bianchini et al., 2022), thereby regulating the water level of the aquifer and the discharge volume to the stream, which is influenced by precipitation intensity, size and distribution of sinkhole. The permeability of sinkholes and conduits typically exhibits multilevel characteristics and varies with scale (Halihana et al., 1999), meaning there are strata structures with different permeabilities, which complicates the flow of water within the karst aquifer and increases the catchment area.

Numerical methods are commonly employed as effective means to accurately simulate karst groundwater movement and assess karst groundwater resources. Shoemaker et al. (2008) proposed a method that discretely embeds conduits, connected by nodes, into the porous media grid (MODFLOW-CFP). This method not only evaluates the water resources of the entire karst aquifer but also considers the geometric shape and distribution of karst conduits on the hydrological processes. Moreover, this methodology has been extensively applied worldwide for estimating karst groundwater flow and water resources (Chang et al., 2015; Qiu et al., 2019; Kavousi et al., 2020; Gao et al., 2020, 2024), as well as in integrated modeling studies coupling SWAT with MODFLOW to investigate groundwater-surface water interactions (Fiorese et al., 2025; Yifru et al., 2024). While MODFLOW-CFP provides robust capabilities for regional-scale karst groundwater simulations, it currently supports only single-phase groundwater flow modeling.

The interaction process between the karst aquifer and stream is also regulated by the dynamic saturation process within the aquifer. The degree of dynamic saturation in different media determines the path and velocity of water flow. Unsaturated aquifers gradually saturate the underlying aquifers under the influence of gravity, while saturated underlying aquifers can cause water to overflow (Worthington, 1991; Huang et al., 2024). In addition, the dynamic saturation processes within the karst aquifer are regulated by factors such as seasonal water level fluctuations, the infiltration and flow of groundwater, and the periodic filling and draining of karst conduits (Huang et al., 2024). It is necessary to couple seepage (porous media) with free flow (conduits and stream) and to describe the dynamic saturation process of the karst aquifer. The Hydrus simulation method based on the Richards equation is capable of simulating variably saturated flow (Dam and Feddes, 2000). However, this approach lacks a built-in conduit flow solution scheme, making it difficult to adequately address the coupling requirements between rapid conduit flow and porous media seepage in karst areas.

Constructing an interaction model between the karst aquifer system and the stream under rainfall event-driven conditions requires coupling free flow and seepage processes while simultaneously supporting two-phase variably saturated flow. (1) The Navier–Stokes (N–S) equations combined with the Darcy equation can effectively couple free flow and seepage processes (Soulaine and Tchelepi, 2016; Carrillo et al., 2020). (2) The Phase Indicator Function for two-phase flow, combined with the phase transition method, can effectively describe the variable saturation process within the karst aquifer (Huang et al., 2024; Zhai et al., 2024). The Darcy–Brinkman–Stokes equations have been utilized to couple seepage flow and free flow (Huang et al., 2024; Nillama et al., 2022; Carrillo et al., 2020). Lu et al. (2023) analyzed a model that integrates fast discharge channels in fractures and conduits with slow seepage in porous media. The results demonstrate that the Darcy–Brinkman–Stokes equations can effectively describe two-phase flow in karst aquifers, and Soulaine (2024) proposed that mixed-scale models based on the Darcy–Brinkman–Stokes equations have strong potential for simulating coupled processes in porous systems.

This study aims to employ a two-phase variably saturated model capable of coupling free flow and seepage flow to reveal the interaction mechanisms between the karst aquifer system and adjacent stream under rainfall infiltration recharge-driven conditions. Specifically, it focuses on further investigating how groundwater saturation variations in different media (e.g., conduits, fractures, matrix) of the karst aquifer system influence inter-media interactions. This research addresses the gap in existing studies where current numerical methods struggle to accurately characterize the collaborative recharge processes among various media within karst aquifer systems. This study employs the Darcy–Brinkman–Stokes equations to model the coupled processes of seepage in porous media and free flow in karst conduit and stream. The Brooks–Corey (BC) and van Genuchten–Mualem (VGM) models are used to characterize the unsaturated seepage in karst media. The Volume of Fluid (VOF) method is applied to monitor the dynamic changes in aquifer saturation. This research elucidates how their saturation dynamics impact the flow exchange among different karst media during precipitation infiltration, and examines the evolving interaction between the karst aquifer and stream under such recharge conditions. Given the complexity of the interaction mechanism between the karst aquifer and stream, this study specifically investigates the impact of four factors on the interaction mechanism: (1) changes in precipitation intensity, (2) different water retention models, (3) multi-stage conduit arrangements, and (4) parameter sensitivity analysis.

2 Materials and methods

The DBS method was employed to couple seepage and free flow, enabling the quantitative characterization of groundwater flow through various media and the interaction processes between the karst aquifer system and adjacent streams.

Unsaturated flow processes within the karst matrix and epikarst zone fundamentally govern the water storage and exchange dynamics. For instance, the shape of the water retention curve determines the amount of water “held” in the matrix at a given suction, thereby controlling the specific moisture capacity and the system's buffer capacity. Meanwhile, the relative permeability function dictates the rate at which hydraulic conductivity decreases as the matrix desaturates. Consequently, these variably saturated processes directly influence the predicted rates of matrix infiltration (during recharge events) and matrix drainage/exfiltration to the conduits (sustaining baseflow), thereby altering the overall storage characteristics and hydrograph response of the karst system.

2.1 Numerical modelling

The numerical model is developed according to the conceptual model of the karst aquifer adjacent to a stream, as depicted in Fig. 1. The model construction incorporates distinct rainfall intensities and temporal rainfall patterns (Fig. 1a and b), while explicitly accounting for characteristic karst geomorphological features including sinkholes, epikarst, and karst conduits. The karst conduit is connected to the epikarst through a sinkhole. The outcrop of the karst spring is located at the end of the karst conduit, directly leading to the stream.

Figure 1 Schematic diagrams of the modelling of the interaction between the karst aquifer (epikarst, sinkhole, karst conduit, PM I (Porous Medium I), PM II (Porous Medium II), and PM III (Porous Medium III)) and stream under dimensionless precipitation intensities (b=3 and b=5). (a) and (a.1) Schematic diagram of the interaction flow between each medium and stream in the early stage of a precipitation event; (b) and (b.1) Schematic diagram of the interaction flow between each medium and stream in the middle stage of a precipitation event. The size of the arrows represents the magnitude of the flow rate, and the direction of the arrows represents the direction of interaction between the two.

Recharge Pathways in a Single Recharge Event: During a single recharge event, precipitation follows two main pathways: a portion directly recharges the adjacent stream, while another portion infiltrates into the epikarst zone (shallow karst system). A fraction of the water stored in the epikarst zone discharges laterally to the stream, while the remaining water disperses vertically through porous media to recharge the deeper porous aquifer. The residual water in the epikarst zone further recharges the karst conduit system via sinkhole point infiltration (Fig. 1a.1).

Conduit Network-Matrix Interaction: Under moderate recharge events, conduits receive water from both sinkhole point recharge and porous media recharge, rapidly transporting it to discharge at karst springs. During intense precipitation events, water in the conduits may temporarily reverse flow to recharge the porous media before returning to the conduits (Bailly-Comte et al., 2010).

Karst Aquifer-Stream Interaction: Lateral recharge from the porous aquifer to the stream requires prior vertical dispersion recharge from the overlying epikarst zone. During a single precipitation event, direct lateral recharge from the epikarst zone and rapid discharge of groundwater from karst springs to the stream cause an earlier stream stage rise. As the stream stage gradually increases, the stream begins to recharge the deeper porous media of the karst aquifer (Fig. 1a.1). Due to the high flow velocity of the stream, its stage declines rapidly, allowing groundwater in the deeper porous media to discharge back into the stream (Fig. 1a.2).

The precipitation influences the dynamic variation process of saturation within porous media, and the water levels in both the karst conduit and the stream experience substantial fluctuations. This variability in water levels is therefore a key driver for the exchange mechanisms between the porous media and the stream. From a hydrological perspective of the watershed, the recharge and discharge processes of karst conduit are controlled by the saturation degree of the surrounding porous media and the water level within the conduit themselves. Based on spatial relationships, the area between the karst conduit and the epikarst is divided into Porous Medium I (PM I) above the conduit, Porous Medium II (PM II) on the sides, and Porous Medium III (PM III) directly below the conduit (Fig. 1a.1). Based on the aforementioned dynamic interaction processes between the karst aquifer system and the adjacent stream, this study constructs the DBS numerical model and employs the CFPv2 (Shoemaker et al., 2008; Giese et al., 2018) to simulate groundwater flow. Through analyses of precipitation intensity variations, multiple precipitation events, different water retention models, multi-level permeability configurations, and parameter sensitivity analyses under repeated rainfall influences, the interaction mechanisms between the karst aquifer system and the stream are elucidated.

2.2 DBS model

2.2.1 Two-Phase Flow Parameter Definition

Assuming that gas and liquid fill the solid pore space, porosity is defined to characterize the percentage of the gas and liquid phases occupying the total pore space.

$$\begin{array}{}\text{(1)}& \mathit{\phi }={\displaystyle \frac{V_{\mathrm{l}}+V_{\mathrm{g}}}{V}}\end{array}$$

In this context, φ represents porosity, V denotes the total volume of the unit [m³], while V_l and V_g correspond to the volumes of the liquid phase (water) and gas phase (air), respectively [m³].

Hirt and Nichols (1981) introduced the VOF method, which employs an additional governing equation to capture fluid motion at free surfaces. Furthermore, the saturation of each phase in the fluid is defined as α_i, where:

• Liquid phase saturation: ${\mathit{\alpha }}_{\mathrm{l}}=\frac{V_{\mathrm{l}}}{V_{\mathrm{g}}+V_{\mathrm{l}}}$.

• Gas phase saturation: ${\mathit{\alpha }}_{\mathrm{g}}=\frac{V_{\mathrm{g}}}{V_{\mathrm{g}}+V_{\mathrm{l}}}$.

Here, the subscripts l and g denote water and air, respectively. Thus, the spatial distribution of water and gas within the porous medium is characterized by porosity φ and phase saturation α_l:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\phi }=\left\{\begin{array}{ll}\mathrm{1}& \text{free\hspace{0.17em} regions}\\ \mathrm{0}<\mathit{\phi }<\mathrm{1}& \text{porous\hspace{0.17em} regions}\\ \mathrm{0}& \text{solid\hspace{0.17em} regions}\end{array}\right.\text{(3)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{l}}=\left\{\begin{array}{ll}\mathrm{1}& \text{water}\\ \mathrm{0}<{\mathit{\alpha }}_{\mathrm{l}}<\mathrm{1}& \text{two-phase\hspace{0.17em} zone}\\ \mathrm{0}& \text{air}\end{array}\right.\end{array}$$

The average fluid density ρ [m³ kg⁻¹] and viscosity μ [m² s⁻¹] within a grid cell are calculated via saturation-weighted averaging:

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\rho }={\mathit{\rho }}_{\mathrm{g}}{\mathit{\alpha }}_{\mathrm{g}}+{\mathit{\rho }}_{\mathrm{l}}{\mathit{\alpha }}_{\mathrm{l}}\text{(5)}& {\displaystyle }\mathit{\mu }={\mathit{\alpha }}_{\mathrm{g}}{\mathit{\mu }}_{\mathrm{g}}+{\mathit{\alpha }}_{\mathrm{l}}{\mathit{\mu }}_{\mathrm{l}}\end{array}$$

where ρ_g is the gas phase density [m³ kg⁻¹] and ρ_l is the liquid phase (water) density [m³ kg⁻¹].

The transport equation for saturation α_i, following Rusche (2002), is expressed as:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial \mathit{\phi }{\mathit{\alpha }}_{\mathrm{l}}}{\partial t}}+\mathrm{\nabla }\cdot \left({\mathit{\alpha }}_{\mathrm{l}}{v}_t\right)+\mathrm{\nabla }\cdot \left(\mathit{\phi }{\mathit{\alpha }}_{\mathrm{l}}{\mathit{\alpha }}_{\mathrm{g}}{v}_{\text{rt}}\right)=\mathrm{0}\end{array}$$

where: v_t is the fluid velocity vector [m s⁻¹], v_rt is the relative velocity between the gas and liquid phases [m s⁻¹].

2.2.2 Governing Equations

To precisely describe groundwater flow through porous media in the karst aquifer system and the free-surface flow processes between conduits and the adjacent stream, this study adopts the DBS equations to characterize immiscible and incompressible two-phase flow in porous media (Nillama et al., 2022; Carrillo et al., 2020; Lu et al., 2023; Huang et al., 2024; Soulaine, 2024). This model provides a unified mathematical framework capable of seamlessly coupling flow phenomena across different scales. This renders it particularly suitable for simulating karst aquifer systems, which are essentially dual-medium systems constituted by both conduits and porous media.

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{v}=\mathrm{0}\text{(8)}& {\displaystyle \frac{\partial \mathit{\phi }{\mathit{\alpha }}_{\mathrm{l}}}{\partial t}}+\mathrm{\nabla }\cdot \left({\mathit{\alpha }}_{\mathrm{l}}\stackrel{\mathrm{‾}}{v}\right)+\mathrm{\nabla }\cdot \left(\mathit{\phi }{\mathit{\alpha }}_{\mathrm{l}}{\mathit{\alpha }}_{\mathrm{g}}{\stackrel{\mathrm{‾}}{v}}_{\mathrm{r}}\right)=\mathrm{0}\text{(9)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathit{\phi }}}& \left({\displaystyle \frac{\partial \mathit{\rho }\stackrel{\mathrm{‾}}{v}}{\partial t}}+\mathrm{\nabla }\cdot \left({\displaystyle \frac{\mathit{\rho }}{\mathit{\phi }}}\stackrel{\mathrm{‾}}{v}\stackrel{\mathrm{‾}}{v}\right)\right)=\\ & -\mathrm{\nabla }\stackrel{\mathrm{‾}}{p}+\mathit{\rho }g+\mathrm{\nabla }\cdot \left({\displaystyle \frac{\mathit{\mu }}{\mathit{\phi }}}\left(\mathrm{\nabla }\stackrel{\mathrm{‾}}{v}+\mathrm{\nabla }{\stackrel{\mathrm{‾}}{v}}^T\right)\right)+F_{\mathrm{c}}+S_{\mathrm{f}}.\end{array}\end{array}$$

Here, t represents the computational time [T], $\stackrel{\mathrm{‾}}{v}$ is the velocity [$L/T$], ${\stackrel{\mathrm{‾}}{v}}_{\mathrm{r}}$ is the relative flow rate of the gas phase to the liquid phase [$L/T$], ρ is the average density of the gas and liquid phases [$M/L^{\mathrm{3}}$], $\stackrel{\mathrm{‾}}{p}$ is the pressure [pa], g is is the acceleration due to gravity (9.81 m s⁻²), μ is the viscosity [$L^{\mathrm{2}}/T$], F_c is the surface tension, and S_f is the resistance source term.

Specifically, within a single set of governing equations, the DBS model is capable of simultaneously describing:

• The high-velocity, free-surface flow within karst conduits;

• The low-velocity seepage flow within the surrounding matrix.

This unification is achieved by strategically incorporating a porosity (φ) and a resistance source term (S_f) into the single momentum conservation equation.

2.2.3 Subdomain Formulation

For the free-flow region and the porous media region, the source terms in the DBS equations adopt distinct forms. Specifically, the source term μk⁻¹ in the two regions can be expressed as (Soulaine, 2024; Huang et al., 2024):

$$\begin{array}{}\text{(10)}& \mathit{\mu }{k}^{-\mathrm{1}}=\mathit{\rho }\mathit{v}{k}^{-\mathrm{1}}+\left\{\begin{array}{ll}\mathrm{0},& \text{free region}\\ k_{\mathrm{0}}^{-\mathrm{1}}{\left(\frac{k_{\text{r,l}}}{{\mathit{\mu }}_{\mathrm{l}}}+\frac{k_{\text{r,g}}}{{\mathit{\mu }}_{\mathrm{g}}}\right)}^{-\mathrm{1}},& \text{porous region}\end{array}\right.\end{array}$$

Here, k₀ is the permeability coefficient determined by the pore structure [m²]. When the permeability is extremely high, this term vanishes, and the DBS equations reduce to the Navier–Stokes (N–S) equations (Eq. 11). Conversely, as permeability decreases, the term $\mathit{\mu }{k}^{-\mathrm{1}}\stackrel{\mathrm{‾}}{v}$ becomes dominant compared to other source terms, causing the DBS equations to asymptotically approach the Darcy equation incorporating gravity and surface tension (Eq. 12).

$$\begin{array}{}\text{(11)}& {\displaystyle }\begin{array}{rl}& \left({\displaystyle \frac{\partial \mathit{\rho }\stackrel{\mathrm{‾}}{v}}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }\stackrel{\mathrm{‾}}{v}\stackrel{\mathrm{‾}}{v}\right)\right)=\\ & -\mathrm{\nabla }\stackrel{\mathrm{‾}}{p}+\mathit{\rho }g+\mathrm{\nabla }\cdot \left(\mathit{\mu }\left(\mathrm{\nabla }\stackrel{\mathrm{‾}}{v}+\mathrm{\nabla }{\stackrel{\mathrm{‾}}{v}}^T\right)\right)+F_{\mathrm{c}},\text{if}\mathit{\phi }=\mathrm{1}.\end{array}\text{(12)}& {\displaystyle }\mathrm{0}=-\mathrm{\nabla }p+\mathit{\rho }g-\mathit{\mu }{k}^{-\mathrm{1}}{v}_{\mathrm{t}}+F_{\mathrm{c}},\text{if}\mathit{\phi }\in ]\mathrm{0},\mathrm{1}[.\end{array}$$

Similarly, the surface tension force F_c and density ρ in the two regions can be expressed as (Huang et al., 2024):

$$\begin{array}{}\text{(13)}& {\displaystyle }{F}_{\mathrm{c}}=\left\{\begin{array}{ll}-\frac{\mathit{\sigma }}{\mathit{\phi }}\mathrm{\nabla }\cdot \left(\frac{\mathrm{\nabla }{\mathit{\alpha }}_{\mathrm{l}}}{\left|\mathrm{\nabla }{\mathit{\alpha }}_{\mathrm{l}}\right|}\right)\mathrm{\nabla }{\mathit{\alpha }}_{\mathrm{l}},& \text{free region}\\ \left[k_{\mathrm{0}}\frac{\left(\frac{k_{\text{r,l}}}{{\mathit{\mu }}_{\mathrm{l}}}{\mathit{\alpha }}_{\mathrm{g}}-\frac{k_{\text{r,g}}}{{\mathit{\mu }}_{\mathrm{g}}}{\mathit{\alpha }}_{\mathrm{l}}\right)}{\frac{k_{\text{r,l}}}{\mathit{\mu }\mathrm{l}}+\frac{k_{\text{r,g}}}{{\mathit{\mu }}_{\mathrm{g}}}}\left(\frac{\partial p_{\mathrm{c}}}{\partial {\mathit{\alpha }}_{\mathrm{l}}}\right)-p_{\mathrm{c}}\right]\mathrm{\nabla }{\mathit{\alpha }}_{\mathrm{l}},& \\ & \text{porous region}\end{array}\right.\text{(14)}& {\displaystyle }\mathit{\rho }=\left\{\begin{array}{ll}{\mathit{\rho }}_{\mathrm{l}}{\mathit{\alpha }}_{\mathrm{l}}+{\mathit{\rho }}_{\mathrm{g}}{\mathit{\alpha }}_{\mathrm{g}},& \text{free regions}\\ k_{\mathrm{0}}\frac{\left({\mathit{\rho }}_{\mathrm{g}}\frac{k_{\text{r,g}}}{{\mathit{\mu }}_{\mathrm{g}}}+{\mathit{\rho }}_{\mathrm{l}}\frac{k_{\text{r,l}}}{{\mathit{\mu }}_{\mathrm{l}}}\right)}{\frac{k_{\text{r,l}}}{{\mathit{\mu }}_{\mathrm{l}}}+\frac{k_{\text{r,g}}}{{\mathit{\mu }}_{\mathrm{g}}}},& \text{porous regions}\end{array}\right.\end{array}$$

Here, σ is the interfacial tension [N m⁻¹], p_c is the capillary pressure [pa], and k_r,g and k_r,l represent the relative permeabilities of the gas phase and liquid phase, respectively.

2.2.4 Relative Permeability Model

Accurate modeling of two-phase flow in porous media is critical in geosciences. Simulating two-phase flow in variably saturated porous media requires precise estimation of the relationship between relative permeability and saturation (Springer et al., 1995).

To characterize the variation in two-phase relative permeability, the effective saturation of the liquid phase must first be defined. This is expressed as:

$$\begin{array}{}\text{(15)}& {\mathit{\alpha }}_{\text{l,e}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{l}}-{\mathit{\alpha }}_{\text{l,r}}}{\mathrm{1}-{\mathit{\alpha }}_{\text{g,r}}-{\mathit{\alpha }}_{\text{l,r}}}}\end{array}$$

where: α_l,e denotes the effective water saturation, α_l and α_l,r represent the water saturation and residual water saturation, respectively, and α_g,r is the residual air saturation.

Relative permeability is a critical parameter in groundwater and related engineering fields (Kuang and Jiao, 2011). The Brooks and Corey (BC) model (Brooks and Corey, 1964) and the van Genuchten model (van Genuchten, 1980) are widely used as representative relative permeability models. The BC model establishes a relationship between relative permeability and effective water saturation as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }{k}_{\text{rg}}={\left(\mathrm{1}-{\mathit{\alpha }}_{\text{l,e}}\right)}^n\text{(17)}& {\displaystyle }{k}_{\text{rl}}={\mathit{\alpha }}_{\text{l,e}}^n\end{array}$$

k_r denotes the relative permeability, where n is a dimensionless coefficient determined by the properties of the porous medium. The Brooks–Corey (BC) model exhibits a sharp discontinuity at the air entry point, which can lead to poor data fitting, particularly for fine-textured soils (Assouline and Or, 2013). The van Genuchten (1980) model addresses this limitation. By incorporating the parameter $m=\mathrm{1}-\mathrm{1}/n$ proposed by Mualem (1976), the modified van Genuchten–Mualem (VGM) model (Parker et al., 1987) is formulated as:

$$\begin{array}{}\text{(18)}& {\displaystyle }{k}_{\text{r,g}}={\left(\mathrm{1}-{\mathit{\alpha }}_{\text{l,e}}\right)}^{\mathrm{0.5}}{\left(\mathrm{1}-{\mathit{\alpha }}_{\text{l,e}}^{\mathrm{1}/m}\right)}^{\mathrm{2}m}\text{(19)}& {\displaystyle }{k}_{\text{r,l}}={\mathit{\alpha }}_{\text{l,e}}^{\mathrm{0.5}}{\left(\mathrm{1}-{\left(\mathrm{1}-{\mathit{\alpha }}_{\text{l,e}}^{\mathrm{1}/m}\right)}^m\right)}^{\mathrm{2}}\end{array}$$

Here, m is a dimensionless parameter.

The selection of permeability equations is critical for appropriate predictions of relative permeability (Yang et al., 2021), indicating that pore tortuosity-connectivity plays a dominant role in groundwater two-phase flow. Therefore, this study conducts simulations and parameter sensitivity analyses for both the Brooks–Corey (BC) and van Genuchten–Mualem (VGM) models.

2.3 CFPv2 model

The CFPv2 model, proposed by Reimann et al. (2014), is an advanced version of MODFLOW-CFP (Shoemaker et al., 2008). It extends functionalities such as flow interactions between conduits and porous media, as well as conduit boundary conditions. CFPv2 integrates with MODFLOW-2005 and employs the following approaches: Laminar Flow in Conduits: Described using the Hagen–Poiseuille equation for discrete conduits within conduit networks. Turbulent Flow: Calculated by combining the Darcy–Weisbach equation with the Colebrook–White equation. Laminar Flow in Fractured Rock Matrix: Simulated via a continuum approach. Detailed technical documentation for MODFLOW-CFP, including groundwater flow simulation methodologies, is provided by Shoemaker et al. (2008). Successful applications and evaluations of the model have been reported in studies such as Gallegos et al. (2013), Reimann et al. (2014), Chang et al. (2019), Gao et al. (2020), and Shirafkan et al. (2023).

2.4 Model Comparison and Numerical Model Construction

2.4.1 DBS Model Conversion and Applicability Assessment

As illustrated in Fig. 2, the Navier–Stokes (N–S) model can resolve fine-scale pore-scale flows and perform high-fidelity simulations. In contrast, the CFPv2 model achieves high computational efficiency and stability by discretizing one-dimensional conduits within porous media. The DBS (Dual-domain Brinkman–Stokes) model combines the advantages of both approaches: By incorporating additional resistance source terms into the N–S equations, it maintains high-fidelity flow resolution in conduits. For porous media, it adopts a Darcy-type flow formulation, significantly reducing computational costs.

Figure 2 Diagram of performance and applicability of different models, (a) N–S model (Navier–Stokes model), (b) DBS model, (c) Schematic diagram of MODFLOW-CFP model solution, (d) Conversion method from DBS equations to N–S equations and Darcy equations.

However, the DBS model operates in three dimensions (3D), requiring grid refinement around conduits and their vicinity to ensure accurate flow resolution. This increases computational load compared to the 1D conduit flow framework of CFPv2. To address this challenge, all simulations in this study were executed on a high-performance server equipped with 64 CPU cores (128 threads) and 256 GB of RAM, which provided the necessary computational power for handling complex 3D meshes.

2.4.2 Model Comparison and Discretization Schemes

To further investigate the effectiveness of the DBS model in addressing interactions between karst groundwater and adjacent streams, this study compares the differences between the MODFLOW-CFP and DBS models. As shown in Fig. 3a.1, the comparison begins with their coupling modes of conduits and porous media from the perspectives of governing equations and grid discretization: MODFLOW-CFP: Groundwater exchange between conduits, porous media, and streams relies on stable hydraulic heads between conduit-porous media and stream-porous media interfaces (Fig. 3a.2). Flow interactions between porous matrix and discrete conduits are linear and driven by head differences (Barenblatt et al., 1960). DBS Model: Groundwater interactions among conduits, streams, and porous media are governed by saturation and pressure gradients between adjacent grid nodes, allowing simultaneous recharge or discharge across interfaces (Fig. 3a.3). However, this requires calculating flux variations across all grids.

Figure 3 (a) Schematic comparison of conduit and porous media coupling modes between MODFLOW-CFPv2 and DBS, (b) DBS model and (c) CFPv2 discretization schemes for karst aquifer systems with riverside models.

Comparison of Stream-Porous Media Interaction Modes: MODFLOW-CFP: Streams are discretized into single grid cells, with exchange fluxes determined by head differences. Fluctuating stream stages are simplified to a uniform water level, and “dry zones” cannot be simulated in porous media (Fig. 3a.4). DBS Model: Media properties (e.g., porosity, permeability) are assigned at grid nodes, and interface values are interpolated. Direct conduit-stream interactions eliminate the need for porous media as an intermediary. Stream geometry can be defined as regular (rectangular) or irregular (Fig. 3a.5). The DBS model employs the Volume of Fluid (VOF) and Front-tracking methods to reconstruct dynamic water-air interfaces, enabling simulation of fluctuating interfaces under sufficiently refined grids.

Discretization Schemes: This study adopts a dynamic programming approach to generate sinkhole and conduit grids, allowing flexible placement of conduits with adjustable diameters and coordinates, enhancing model adaptability (contrasting fixed conduit positioning in studies like Kavousi et al., 2020; Pardo-Igúzquiza et al., 2018; Li et al., 2023).

DBS Discretization (Fig. 3b): The epikarst layer thickness and stream location are defined. Regions are divided into free-flow zones (streams, sinkholes, conduits) and porous media. Free-flow zones use locally refined grids to capture micro-scale variations in water levels and interfaces. Porous media zones adopt gradually coarsening grids (edge cells twice the size of conduit-adjacent cells), balancing accuracy and computational efficiency. Permeability is graded, decreasing outward from conduits to reflect dissolution effects.

CFPv2 Discretization (Fig. 3c): Conduits are embedded in porous media and directly connected to streams. Domain dimensions: $\mathrm{200}\mathrm{m}\times \mathrm{200}\mathrm{m}\times \mathrm{30}\mathrm{m}$ (length × width × thickness). Groundwater flows from porous media to conduits and discharges into streams (Fig. 11a.1).

Porous media: Homogeneous, initial head =10 m, no-flow boundaries. Conduits: Diameter =1 m, roughness =0.01 m, wall interaction parameter =25 m s⁻¹, outlet collocated with stream grid. Initial conditions: Spring discharge =0, conduit node elevation =1 m, water temperature =20 °C. Boundary conditions: Rainfall recharge at the top, total simulation time =45 000 s, MODFLOW-CFP stress periods =1 min.

2.5 Rainfall Infiltration Recharge Boundary

The upper boundaries of both the DBS and CFPv2 models are defined as transient natural precipitation boundary conditions. In this study, the rainfall infiltration recharge boundary condition is formulated as follows (Huang et al., 2024; Chang et al., 2015):

$$\begin{array}{}\text{(20)}& I\left(t\right)={\displaystyle \frac{b}{\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}}\sum e^{-\frac{{\left(\frac{t_i-\mathit{\mu }}{a}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\end{array}$$

Here, t_i denotes the time of the ith rainfall event, and I(t) represents the rainfall intensity at time t. According to Chang et al. (2015), the parameters μ, σ², and a are set as constants (90, 1.5, and 20, respectively). Variations in rainfall intensity during the infiltration recharge process, along with the total amount and peak intensity of the event, are controlled by adjusting the dimensionless scaling parameter b.

3 Results

3.1 Interaction process between the karst aquifer and stream under precipitation infiltration recharge

3.1.1 Karst aquifer and stream interactions under varying precipitation intensities

The changes in hydrological process curves, water level fluctuations, and their differences during the interaction between karst media and stream under different precipitation intensities are shown in Fig. 4. In the early stage of precipitation, the flow in the stream primarily originates from direct precipitation recharge and lateral groundwater recharge from epikarst (Fig. 4a). As the water level in the stream gradually rises, the flow not only continues downstream but also begins to recharge the karst aquifer, particularly the PM II. The peak recharge to PM II coincides with the peaks of epikarst recharge to the stream (Epikarst in Fig. 4) and direct precipitation recharge (P-River in Fig. 4). Therefore, the interaction process between the karst aquifer and stream during the early precipitation stage is significantly influenced by lateral groundwater discharge from the epikarst and the direct precipitation recharge. As groundwater recharge from epikarst to the stream declines (Fig. 4a), groundwater moves downward through the epikarst to PM I, and begins to gradually recharge the stream. Due to the low permeability of the epikarst, lateral discharge from PM I to the stream will be delayed. During this process, the discharge volume of PM I exhibits two distinct peaks. The first peak is due to the recharge of groundwater from the epikarst, while the second peak is caused by the gradual saturation of PM II and the karst conduit, with a proportion of groundwater overflowing from PM I and discharging laterally to the stream. After the end of precipitation recharge, the hydrological process curve of PM I rapidly declined, and the discharge volume of the karst conduit, PM III and PM II gradually increase, causing the water level in the stream to rise (Fig. 4d). When the water level in the stream gradually exceeds that of PM I, the stream begins to gradually recharge PM I. The karst conduit, PM II and PM III continue to discharge to the stream during this stage due to higher internal water pressure, forming a local hydrological cycle with the upper layer. In the late stage of precipitation, the hydrological process of the stream primarily shows a gradual decline in baseflow.

Figure 4 Hydrological process curves of each medium in the karst aquifer and stream for different precipitation intensities: (a) b=3, (b) b=5, (c) b=7. Water level changes and differences in water levels in the karst aquifer and stream for different precipitation intensities: (d) b=3, (e) b=5, (f) b=7.

As depicted in Fig. 3b and c, the recharge and discharge dynamics between the karst aquifer and stream across different media shift notably with escalating precipitation intensity. The recharge volumes from the stream to PM I and PM II both decrease. The reduction in the recharge to PM II from the stream is primarily due to the acceleration of groundwater movement downward as precipitation intensity increases, causing groundwater to move more rapidly to the bottom of the karst aquifer, thereby recharging PM II. Consequently, part of pore space that should have been recharged by the stream is instead recharged from PM I downward. The decrease in the recharge to PM I can be attributed to its high internal saturation level and the rise in water level. On the other hand, the water level in the stream does not significantly exceed that of the upper aquifer, making it difficult for the stream to effectively recharge the aquifer. Due to the reduced recharge volume to the aquifer, the discharge from the stream is partially lower than the discharge from the epikarst during the early stage of the hydrological process.

With changes in precipitation intensity (b=3, 5, and 7), the water level variations and their differences between the karst aquifer and stream exhibit complex dynamic characteristics (Fig. 3d–f). During the early stage of precipitation, despite the increasing water level difference, the discharge from the stream to the aquifer is gradually decreasing (as shown by the negative values for PM I and PM II in Figs. 4a, 3b and c). This phenomenon indicates that water level is not the only factor controlling the interaction between the karst aquifer and stream; changes in the degree of saturation also play a significant role. As shown in Fig. 4d, under low precipitation intensity, the water level difference between the karst aquifer and stream is often greater than the water level of the stream during the middle and later stages of precipitation. However, as precipitation intensity increases, the water level difference tends to decrease (Figs. 4b and 3c). This change is primarily due to the increased precipitation intensity leading to a faster saturation of the karst aquifer, thereby limiting the ability of the stream to recharge the aquifer. After the middle stage of precipitation, the interaction between the stream and the upper part of the aquifer gradually intensifies, while the lower part of the aquifer discharges to the stream (Fig. 4a). Due to the gradual decrease in water level difference, it is difficult for the stream to effectively recharge the aquifer. In this process, the interaction between the aquifer and stream is controlled by the dynamic changes in saturation.

Based on the comparison between DBS and Modflow-CFPv2 results in Fig. 4a–c, the CFPv2 model exhibits a single-peak hydrograph with exponential recession characteristics, failing to capture flow process line disturbances caused by multi-media interactions. Under precipitation intensities b=3 and 5, the CFPv2 model shows an immediate rapid increase in stream discharge during early stages rather than gradual enhancement, though total discharge and baseflow during later stages remain comparable (as shown in Table 3). Specifically, for b=3, the peak stream discharge in Modflow-CFPv2 occurs at 2520 s, earlier than in the DBS model. This discrepancy arises because the precipitation recharge package in CFPv2 directly elevates water levels, whereas the DBS model simulates a gradual vertical infiltration process along the z axis. Lower precipitation intensity reduces groundwater infiltration rates and prolongs water table replenishment time, consequently delaying lateral discharge timing. At b=7, both models exhibit comparable first discharge peaks, but the DBS model generates a secondary peak through overflow effects that rapidly recedes after overflow cessation. In contrast, CFPv2 demonstrates smooth exponential recession without secondary features due to its simplified vertical stratification that neglects multi-component interactions.

Table 1 Variable Definition Table.

Table 2 Different parameter used in Models.

Table 3 Comparing DBS and MODFLOW results for key variables.

The comparable results between DBS and Modflow-CFPv2 models under variable recharge conditions demonstrate the reliability and stability of DBS in simulating karst aquifer systems. Although the DBS model captures more interaction details, it requires greater computational resources. The absence of overflow mechanisms and multi-media interactions in CFPv2 leads to simplified discharge recession patterns that fail to reflect intense component interactions within the system. This comparative analysis highlights the DBS model's advantages in characterizing complex conduit-stream-aquifer interactions while acknowledging its computational demands.

It is self-evident that changes in precipitation intensity significantly affect the recharge and discharge processes between the karst aquifer and stream. The water levels and saturation degrees of the respective media act as core controlling factors that jointly influence the interactive dynamics between the aquifer and stream. To gain a deeper understanding of these influencing factors and their interaction mechanisms, and to further elucidate the interaction process mechanisms between the karst aquifer and stream, this study focuses on the hydrological interaction process between the two during the early stage of precipitation.

3.1.2 Interaction process between the karst aquifer and stream during early stage of precipitation

Figure 5 illustrates how the interaction volume between the epikarst, porous media, and stream varies under different precipitation intensities. As shown in Fig. 5a, at a precipitation intensity b=3, the contribution ratios of the epikarst, PM I, and PM II to the recharge of the stream are similar. This indicates that during the early stage of precipitation, the recharge effects of each medium on the stream are relatively balanced. Since groundwater vertically recharges the underlying aquifer through the epikarst, the discharge peak of PM II is relatively delayed compared to the epikarst and PM I.

Figure 5 Interaction process of epikarst, porous media, and stream for different precipitation intensities: (a) b=3, (b) b=5, (c) b=7.

As the precipitation intensity increases (b=5), the contribution ratios of the epikarst, PM I, and PM II to the recharge of stream experience significant changes (Fig. 5b). Upon comparing Figs. 5a and 4b, it is evident that an increase in precipitation intensity leads to higher discharge volumes for both PM I and PM II, with PM II experiencing a more pronounced rise. Additionally, the peaks of their discharges occur earlier. The first peak of PM I is primarily caused by infiltration recharge from precipitation. With the increase in precipitation intensity, the infiltration velocity accelerates and the recharge volume increases, leading to a larger discharge volume and an earlier peak for PM I (vertical recharge peak). Groundwater continues to move downward from PM I, and the saturation of PM II rises, allowing more groundwater to overflow and discharge through PM I, thereby generating the second peak (overflow peak). For PM II, as discussed in Sect. 3.1, increase in saturation reduces the recharge from stream, but the discharge volume increases gradually after the middle stage of precipitation, and its contribution to the recharge of the stream becomes dominant among the three. This is due to the increased precipitation intensity, which allows PM II to receive more vertical recharge, enhancing its discharge capacity. When the precipitation intensity continues to increase (b=7, Fig. 5c), PM II gradually reaches saturation. According to the analyses in Sect. 3.1, the ability of PM II to receive recharge is limited by its own saturation level, making it difficult to receive vertical recharge. Therefore, despite the increased precipitation intensity, the discharge volume of PM II does not increase significantly. Conversely, due to the influence of the saturation state of the underlying aquifer medium, the second peak (overflow peak) of PM I is more pronounced, indicating a more evident overflow phenomenon. Under higher precipitation intensity, the recharge contribution of PM I to the stream dominates.

Thus, variations in precipitation intensity notably influence the interaction volume between the karst media and stream. As precipitation intensity increases, the discharge volume and peak values of each medium are altered. Specifically, the two peaks of PM I show sequential changes in intensity, which are modulated by the saturation levels of the adjacent media.

3.1.3 Dynamic interaction processes between various media within a karst aquifer system

The DBS model, leveraging its fine grid resolution and two-phase flow simulation capability, can accurately capture the interactive processes between various media (e.g., saturated-unsaturated zones, conduit-stream systems) influenced by dynamic saturation processes during precipitation infiltration recharge. As the interactions between adjacent media are governed by variations in saturation levels, the numerical results under rainfall intensity b=5 are selected for further analysis of dynamic inter-media interactions.

As shown in Fig. 6, the Darcy–Brinkman–Stokes model clearly demonstrates the changes in the saturation levels of epikarst, porous media, and the karst spring; the saturation fields and the interaction between various media at 4000, 6105, and 7363 s; the interaction amounts between epikarst, porous media I, II, III, and the stream. From Fig. 6a.1, it can be seen that the saturation level of epikarst rises and declines earliest, but the saturation level is relatively low, and it is in a completely unsaturated flow state. Porous media I and III rise synchronously before 5000 s, while porous media II and the karst spring rise rapidly at 4611 s. At 7409 s, the karst spring and porous media I successively enter the decline stage. Due to the rapid drainage of the conduit, the saturation level decreases. The saturation level of the karst spring decreases faster than that of porous media I and intersects with porous media I at 9670 s.

Figure 6 For the Darcy–Brinkman–Stokes model: (a.1) Variations in the saturation levels of epikarst, various porous media, and the karst spring. (a.2) Saturation fields and the interaction among different media at 4000, 6105, and 7363 s. (a.3) Interaction volumes between epikarst, porous media I, II, and the stream. (a.4) Interaction volumes among the karst spring, porous media III, and the stream.

Combining Fig. 6a.2 with other sub-figures, the stages with obvious interactions among porous media can be divided into the infiltration stage (green), the overflow stage (red), and the recession stage (blue). During the infiltration stage from 4000 to 4611 s, as shown in Fig. 6a.2.1, epikarst vertically replenishes PM I and infiltrates downward. However, the infiltrating water does not reach the lower media. Meanwhile, the saturation levels of porous media II, III, and the conduit gradually increase (see Fig. 6a.1). Combining with Fig. 6a.3, it can be seen that epikarst laterally replenishes the stream, and quickly drops to the bottom of the riverbed due to gravity. At this time, the lower aquifer system (porous media II, III, and the conduit) is in a dry state, so the stream replenishes the lower aquifer. The amount of recharge received by PM III and the conduit is less than that of PM II (analyzed by combining Fig. 6a.3 and a.4), but their saturation levels increase faster. There are two reasons for this situation: First, the bottom elevation of the conduit is 1 m, and the water level of the stream needs to submerge the 1 m water level before it can recharge the conduit. Second, PM III is not only replenished by the stream, but also the sinkhole diverts the groundwater in epikarst and PM I to the conduit (the sinkhole flow velocity and saturation as shown in Fig. 6a.2.1), and then replenishes PM III. As the lower aquifer media gradually tends to be saturated with rainfall recharge, as shown in Fig. 6a.2.2, porous media II and III tend to be saturated (see Fig. 6a.2.1). Due to the weak compressibility of water, after the upper part infiltrates and replenishes PM I, it tends to laterally replenish the stream from the interface between PM II and stream. As the saturation level of PM I gets higher, the lateral recharge to the stream becomes more significant, showing an obvious overflow state. The depression between the two peaks is caused by the rapid rise of the stream water level. During the flood peak stage, the discharge from porous media to stream decreases. At the same time, the rise of the stream water level makes it difficult for the lower porous media to replenish the stream, and PM II tends to be saturated, making it difficult to replenish PM I. During this stage, the flow between porous media I and II is in a dynamic equilibrium state. As shown in Fig. 6a.2.3, during the recession stage, the rainfall infiltration intensity decreases rapidly. Under the action of gravity, the groundwater vertically replenishes PM I, the conduit, and PM II successively recedes. And the water level of the stream drops rapidly (see Fig. 3e). The groundwater tends to be discharged to the stream through PM I and the karst spring. PM I is replenished by PM II on the one hand and discharges to the stream on the other hand. Therefore, during a single rainfall event, during the infiltration stage, part of the amount of water replenished from epikarst to the stream is discharged, and other part is redirected to replenish the lower porous media; during the overflow stage, the stream is mainly replenished through the karst conduit and PM II. PM I and the stream are in a dynamic equilibrium state. During the recession stage, the porous media act as the main medium to replenish the stream.

As shown in Fig. 6a.4, the karst spring reaches its peak at 7409 s. This is due to the rainfall infiltration, the recharge from PM I, and the subsequent discharge to the stream. As the storage volume decreases, the amount of recharge from the karst spring to the stream decreases. A trough appears at 11 642 s. This is because as the water level of the stream drops, groundwater is more easily discharged into the stream. However, as the overall storage volume continues to decline, after a peak appears at 13 057 s, it enters a complete recession stage. Affected by the decline of the stream water level, the discharge from PM III to the stream gradually increases during the recession stage. Combining with Fig. 6a.1, it can be seen that while PM III is discharging, its saturation remains at level I continuously, indicating that the conduit continuously supplies water vertically to PM III.

Under the recharge of rainfall infiltration, the interaction process between the karst aquifer affected by epikarst, sinkholes, conduit and the stream shows dynamic changes in terms of staged characteristics, main interaction media, and the dynamic equilibrium process among different media. The accurate simulation of the above complex processes depends on the support of a 3D two-phase DBS model

3.2 Impact of multiple precipitation events on the interaction process between the karst aquifer and stream

Rainy seasons typically experience multiple precipitation events, during which differences in precipitation peaks, durations, and cumulative precipitation events can all impact the interaction process between the karst aquifer and stream. Does the groundwater stored in the porous media of the karst aquifer system during the initial rainfall event influence the interactions between multi-component media during subsequent precipitation episodes?

Based on understanding the interaction mechanism of a single precipitation event, this study further analyzes the impact of multiple precipitation events on the interaction process. Figure 7 shows the changes in water level of stream under continuous precipitation events. When the intensities of two consecutive precipitation events remain constant, the water level of stream reaches both the highest and the lowest points, indicating that the water level is related to the total precipitation intensity. Even with different intensities of the first precipitation event (b₁=3 and b₁=5), the trend of the water level changes in stream is consistent (Fig. 7① and ④). After the first precipitation event, the karst aquifer receives infiltration recharge from the precipitation and can store part of the water, so the water level of stream will be higher during the second precipitation event, and the greater the intensity of the second precipitation event, the higher the water level of stream (Fig. 7① and ②, or ③ and ④). This indicates that the intensity of the second precipitation event determines the amount of recharge from each medium to stream. Therefore, when the intensity of the first precipitation event is the same, the amplitude of the water level change in stream during the second precipitation event is only related to the intensity of the second precipitation event. When the intensity of the second precipitation event is the same, the storage capacity of the karst aquifer during the first precipitation event determines the amplitude of the water level change in stream during the second precipitation event. When the total precipitation intensity is the same (Fig. 7② and ③), if the intensity of the first precipitation event is lower than that of the second one, the amplitude of the water level change in stream is higher, and vice versa. This is because, in the case of two consecutive precipitation events, part of the precipitation infiltrates and recharge the storage during the first event, and the other part is discharged to stream through the aquifer. Combining Fig. 4d and e, during the first precipitation event, the water level in the porous medium rises and stores a proportion of water, but the discharge volume to stream is greater when the precipitation intensity is higher (b₁=5) compared to when it is lower (b₁=3, Fig. 4a and b). When the second precipitation event occurs, due to the similar saturation levels of the karst aquifer, the greater the intensity of the second precipitation event, the larger the amount of groundwater recharged to stream through the aquifer, and the more pronounced the amplitude of the water level in stream.

Figure 7 Water levels in stream for two consecutive precipitation events with first and second precipitation intensities: Scenario ① b₁=3 and b₂=3; Scenario ② b₁=3 and b₂=5; Scenario ③ b₁=5 and b₂=3; and Scenario ④ b₁=5 and b₂=5, respectively.

Figure 8 illustrates the hydrological process curves of the stream during two consecutive precipitation events, as well as the interaction processes between the various media of the karst aquifer and stream. Under different precipitation intensities, the various media of the karst aquifer recharge the stream with varying intensities, resulting in significant fluctuations in the water level of stream. Based on Figs. 8a and 7② and ④), it can be observed that under two consecutive precipitation events, when the intensity of the second precipitation event is equal to or greater than the first, the stream hydrograph exhibits more pronounced fluctuations. The comparison between the DBS model and MODFLOW-CFPv2 model under different b1 parameter combinations demonstrates distinct characteristics in streamflow hydrographs: the DBS model shows higher peak discharge with greater fluctuations, while the MODFLOW-CFPv2 model displays relatively smoother discharge variations. Notably, under the second precipitation event, the MODFLOW-CFPv2 model exhibits delayed peak elevation timing. Furthermore, its recession phase still follows an exponential decay pattern, failing to capture the rapid interactive response between multi-media systems during successive precipitation events. As shown in Fig. 8b, the epikarst discharges quickly and is not easily affected by multiple precipitation events. However, when the intensity of the first precipitation is high and the intensity of the second precipitation is the same (① and ③), the discharge volume of the epikarst to stream during the second precipitation period is slightly larger. When the intensity of the first precipitation is different and the intensity of the second precipitation is the same (Fig. 8c② and ④), the discharge volume of groundwater through karst conduit to stream during the second precipitation period is almost the same. This is because karst conduit discharge quickly, and the storage volume of the conduit during the first precipitation period has little impact on the storage volume during the second precipitation period. Therefore, combining with Fig. 7, it is known that the storage effect of the karst aquifer mainly occurs in the porous medium, and it also indicates that relying solely on changes in the water level of stream makes it difficult to clearly determine the storage volume of the porous medium and conduit during the first precipitation event, and their respective impacts on the second precipitation period (Fig. 7). When the intensity of the second precipitation is higher (Fig. 8c②, ③ and ④), the discharge volume of the porous medium (PM II) to stream does not increase significantly. This is because the intensity of the second precipitation is larger, causing the water level of stream to rise (Fig. 7), making it difficult for the porous medium (PM II) to recharge stream.

Figure 8 (a) Hydrological process curves of the stream; (b) Discharge process of groundwater through the epikarst to the stream; (c) Discharge process of groundwater through the karst conduit to the stream; (d) Discharge process of porous media (PM II) to the stream, for two consecutive precipitation events with first and second precipitation intensities: Scenario ① b₁=3 and b₂=3; Scenario ② b₁=3 and b₂=5; Scenario ③ b₁=5 and b₂=3; and Scenario ④ b₁=5 and b₂=5, respectively.

Therefore, under the influence of two consecutive precipitation events, the greater the total precipitation intensity, the larger the discharge volume of the karst aquifer to stream. The storage effect of the karst aquifer occurs in the porous medium and affects subsequent precipitation processes. The lower-level porous medium (PM II), due to the high water level and large fluctuations of stream, is more difficult to recharge stream, and the recharge from stream mostly comes from overflow supply from the media in other layers.

3.3 Effects of Water Retention Characteristics on Karst Aquifer-Stream Interactions

The external recharge of the system significantly influences the interaction processes among different media. This study further investigates how the inherent hydrogeological properties of karst systems affect these interactive processes. Variable saturated flow in the karst vadose zone plays a critical role (Dvory et al., 2018), where the water retention characteristics of porous media govern unsaturated flow dynamics. However, the CFPv2 model struggles to simulate variable saturation processes. This paper compares the DBS model results with two distinct experimental datasets to elucidate the advantages and limitations of the DBS approach in simulating variable saturated flow. This study selected the experiments by Warrick et al. (1985) and Vauclin et al. (1979) because, although these physical experiments have fewer data points (compared to modern numerical simulations), they clearly demonstrate the transient evolution of pressure head or water table position. This is both necessary and sufficient to validate the DBS model's capability in handling variably saturated flow – a capability that CFPv2 lacks.

Case 1. A typical unsaturated-unsteady seepage problem in sandy clay loam (Warrick et al., 1985), where the soil hydraulic properties are provided by the international UNSODA database (Leij et al., 1996). Key parameters include: $k=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, α_s=0.363, α_r=0.186, and n=1.53. The model consists of a vertical soil column (1 m thickness) with an initial pressure head of −8 m across the domain. The top boundary is set to a pressure head of 0 m to simulate free surface infiltration.

Case 2. A 2D laboratory infiltration experiment by Vauclin et al. (1979), widely used for evaluating saturated-unsaturated unsteady seepage models. The soil slab measures 2.00 m in height, 6.00 m in width, and 0.05 m in thickness, with an impermeable base and free drainage boundaries on both sides. Initially, the water table is set at 0.65 m. A central 1.00 m section of the top boundary receives uniform precipitation at 0.148 m h⁻¹ for 8 h, during which free surface evolution is monitored. Soil hydraulic properties are described using the van Genuchten–Mualem model with parameters: k=0.35 m h⁻¹, α_s=0.30, α_r=0.01. Due to symmetry, the DBS model simulates the right half of the domain.

As shown in Fig. 9, the DBS model demonstrates strong agreement with both experimental datasets, highlighting its capability to capture spatiotemporal variations in water-air two-phase flow. Comparative analysis between DBS simulations and experimental data not only validates model reliability but also enhances understanding of soil moisture transport mechanisms. This provides critical support for simulating interactions between karst aquifers and adjacent streams.

Figure 9 Comparison between the DBS model and experimental results from (a) Warrick et al. (1985) and (b) Vauclin et al. (1979).

Based on the well-validated two-phase flow DBS model, this study analyzes the impacts of different water retention models on interactive flow between media. Figure 10 presents the hydrograph curves under different water retention model parameters (BCn=3, 2.5, 2 and VGMm=0.85, 0.8) for (a) stream, (b) karst spring, (c) epikarst, (d) PM I, (e) PM II, and (f) PM III. Figure 10c.1 illustrates the parameter effects on porous media morphology, where n≥2 and higher n values indicate more heterogeneous pore space and complex structures. Figure 10d.1 compares water retention curves between BC and VGM models.

Figure 10 Hydrological process curves under different water retention model parameters (BCn=3, 2.5, 2 and VGMm=0.85, 0.8) for (a) stream, (b) karst spring, (c) epikarst, (d) PM I, (e) PM II, and (f) PM III. Subplots (c.1) and (d.1) show the schematic diagram of parameter effects on porous media morphology and the water retention curves of the BC and VGM models, respectively.

Combining Fig. 10a and b, in the BC model, increasing n values progressively reduce hydrograph curves of stream and karst spring, attributed to irregular pore media impeding groundwater flow and reducing discharge. In the VGM model, decreasing m values (equivalent to increasing n) enhance pore structure irregularity, similarly lowering hydrograph curves. As shown in Fig. 10c, epikarst discharge increases with higher n values due to its low permeability (K0) during relative permeability correction, facilitating enhanced groundwater discharge through epikarst to the stream.

From Fig. 10d and e, larger n values correspond to decreased epikarst-stream discharge and increased downward recharge to porous media, thereby enhancing stream recharge from PM I and II. Integrating Fig. 10c and e, reduced epikarst-stream hydrographs with higher n values lead to diminished stream-porous media recharge. Figure 10f demonstrates that PM III is primarily influenced by conduit flow and shows minimal sensitivity to n and m parameters.

Figure 10d.1 displays saturation variations derived from two karst groundwater retention models: Brooks–Corey (BC) model (Eqs. 16 and 17) and van Genuchten–Mualem (VGM) model (Eqs. 18 and 19). For identical infiltration periods, BC model predicts higher moisture retention than VGM. The BC model emphasizes static water retention in karst media, while VGM prioritizes dynamic groundwater transport and distribution. The VGM model predicts longer groundwater migration distances, suggesting greater sensitivity in simulating karst groundwater diffusion and infiltration processes. These differences hold significance for unsaturated two-phase flow dynamics and accurate prediction of groundwater migration paths in karst aquifer systems.

Furthermore, discrepancies exist between BC and VGM models in simulating saturation variations (Fig. 10d.1), manifesting as distinct saturation degrees and groundwater migration distances under identical conditions. Therefore, selecting appropriate models based on lithological characteristics is crucial for precise description and prediction of two-phase flow in karst groundwater systems.

3.4 Impact of multi-stage permeability and porosity arrangement on the interaction process between the karst aquifer and stream

In this study, the “multi-level conduit” configuration is our model's conceptualization of the “nested hydraulic discontinuities” (Halihana et al., 1999) inherent to karst, representing the spectrum of heterogeneity created by the co-existing matrix, fracture, and conduit flow components. By comparing the multi-level and single-level conduit configurations, the results show that the configuration choice did not induce significant changes in the hydrological processes of the epikarst (Fig. 11c), PM I (Fig. 11d), and PM II (Fig. 11e). In these media, the “M” and “S” hydrographs are nearly identical. However, the impact of the multi-level configuration was significant for the main stream (Fig. 11a), the total karst system discharge (Fig. 11b), and PM III (Fig. 11f). In all these cases, the multi-level (M) configuration resulted in a visibly higher and earlier peak discharge compared to the single-level (S) configuration. As shown in Fig. 11a, when multi-level conduit arrangements are adopted, the peak of stream hydrological process increases, indicating that multi-level conduit arrangements enhance the recharge volume of stream. However, during the recession phase, the flow under multi-level conduit arrangements is relatively low. This is because multi-level conduit collects a proportion of the flow that should have been contributed by the later stage matrix recession and discharge it to stream, thereby affecting the peak of the recession process. As shown in Fig. 11b, under multi-level conduit arrangements, sinkhole can absorb more water and discharge it through karst conduit. This indicates that multi-level conduit arrangements can more effectively play their roles in water absorption and discharge during heavy precipitation events. However, in the case of lower precipitation intensity in the early stage, the water absorption priority of multi-level conduit is not fully manifested. By comparing Fig. 11c–e, it is found that multi-level conduit arrangements have no significant impact on the hydrological processes of the epikarst and porous media (PM I and PM II). This suggests that multi-level conduit arrangements mainly affect the interaction between the karst conduit and stream, with relatively little impact on other media. The hydrological responses of the karst conduit and PM II under multi-level conduit arrangements are shown in Fig. 11f and b. Under multi-level conduit arrangements, the discharge volume of the karst conduit significantly increases. At the same time, due to the increase in karst conduit flow, PM II also receives more recharge, leading to a corresponding increase in the discharge volume of this portion of porous media to stream. This further indicates that multi-level conduit configurations can notably influence the hydrological processes of stream and karst conduit under specific precipitation intensities, with minimal effects on other media.

Figure 11 Impacts of single-stage and multi-stage conduit hydrological process changes in various media of the karst aquifer for a precipitation intensity b=5.

4 Uncertainty Analysis and Discussion

The multi-level conduit configuration inherently affects multi-media interactions by simultaneously altering permeability, conduit diameter, and porosity parameters. This study will further conduct sensitivity analyses on individual variables to investigate their impacts on the vulnerability of karst aquifer systems.

4.1 Impacts of Conduit Diameter and Geometry on Interactions Between Karst Aquifer Systems and Streams

Figure 12 presents hydrographs under conditions of circular conduits with varying radii (r=0.2, 0.3, 0.3, and 0.5 m) and square-section conduits (r=0.5 m) for (a) stream-connected flow, (b) karst spring discharge, (c) epikarst flow, (d) PM I (PM I), (e) PM II, and (f) PM III. Figure 12c.1 illustrates different conduit cross-sectional shapes to analyze their impacts on the interactive flow between karst aquifer systems and adjacent streams.

Figure 12 Hydrological process curves for (a) stream, (b) karst spring, (c) epikarst, (d) PM I, (e) PM II, and (f) PM III under conditions of circular conduits with radii r_c=0.2, 0.3, 0.3, and 0.5, and square-cross-section conduits with $S-r_{\mathrm{c}}=\mathrm{0.5}$. Subplot (c.1) shows a schematic diagram of different conduit cross-sectional shapes.

As shown in Fig. 12a, larger conduit radii correspond to higher initial discharge peaks and shorter peak arrival times, indicating enhanced porous medium recharge and faster fluid transmission through larger conduits. Notably, the square-section conduit ($s-r_{\mathrm{c}}=\mathrm{0.5}$) exhibits higher peak discharge than its circular counterpart (r_c=0.5) due to its surplus cross-sectional area accommodating greater fluid discharge under identical nominal radii.

Figure 12b demonstrates that karst spring peak discharge increases with conduit radius. At r=0.5 m, the square-section conduit ($s-r_{\mathrm{c}}=\mathrm{0.5}$) achieves higher peak discharge than the circular conduit (r_c=0.5), but displays lower recession flow. This occurs because identical precipitation infiltration recharge leads to greater porous medium storage depletion during peak periods in square conduits, subsequently reducing porous medium-to-conduit recharge during baseflow recession.

Combined analysis of Fig. 12c–e reveals that conduit radius variations do not significantly affect epikarst hydrographs or PM I/II hydrographs. However, square-section sinkholes modify flow patterns: epikarst hydrographs show lower values under square conduits, while PM I/II hydrographs exhibit higher values due to enhanced epikarst groundwater collection in square cross-sections, increasing recharge to PM I/II.

Figure 12e indicates that larger conduit radii correspond to lower negative values. Combined with Fig. 12a, this demonstrates that increased stream recharge through larger conduits elevates both stream peak discharge and water levels, thereby enhancing porous medium-stream interactions. Similarly, Fig. 12f shows that larger conduit radii increase karst spring discharge and PM III hydrograph elevation through enhanced gravity-driven groundwater recharge.

Conduit geometry (radius and shape) constitutes a critical factor in karst aquifer hydrological modeling. Larger circular conduits accelerate peak discharge arrival and amplify stream-connected flow peaks and karst spring discharge. Square-section conduits outperform circular equivalents in peak discharge capacity under identical nominal radii due to cross-sectional area advantages. Enlarged conduits intensify porous medium-stream interactions and amplify PM III recharge through gravitational effects. Comprehensive consideration of conduit geometry impacts on hydrological elements is essential for improving model accuracy and reliability in simulating karst aquifer-stream interaction processes.

4.2 Influence of Permeability on the Interaction Processes Between Karst Aquifer Systems and Streams

The permeability of the epikarst directly controls the ease of fluid infiltration from the surface into the conduit system. Figure 13 illustrates the hydrological process curves under different epikarst permeability coefficients ($K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{6}},{\mathrm{10}}^{-\mathrm{7}},{\mathrm{10}}^{-\mathrm{8}},{\mathrm{10}}^{-\mathrm{9}}$; when $K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{9}}$, the permeability matches that of porous media, rendering the epikarst incapable of rapid groundwater leakage) for: (a) stream, (b) karst spring, (c) epikarst, (d) PM I, (e) PM II, and (f) PM III. This aims to reveal how epikarst permeability regulates groundwater flow patterns in complex conduit systems and intermedia interactions.

Figure 13 Hydrographs under different epikarst permeability conditions ($K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{6}}$, $K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{7}}$, $K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{8}}$, $K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{9}}$) for: (a) stream, (b) karst spring, (c) epikarst, (d) PM I, (e) PM II, (f) PM III. Subfigure (c.1) shows a schematic diagram of media interactions under varying epikarst permeability conditions.

As shown in Fig. 13a, under high epikarst permeability ($K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{6}}$): the discharge curve rises rapidly to a peak of ∼4.5 m³ s⁻¹ followed by a sharp decline. This indicates that high permeability enables rapid groundwater leakage from the epikarst to the stream, causing swift flow increases. Peak stream discharge diminishes with decreasing permeability. High permeability reduces flow resistance, facilitating faster fluid entry into the conduit system and generating sharp discharge peaks, while low permeability increases resistance, resulting in gradual fluid release and broader, lower discharge curves.

Figure 13b demonstrates that epikarst permeability differences from porous media have minimal impact on conduit flow. However, when epikarst permeability equals that of porous media ($K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{9}}$), the peak discharge at the karst spring decreases while maintaining identical baseflow recession characteristics. Combining Fig. 13c and c.1, higher epikarst permeability enhances lateral discharge to the stream. At $K_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{9}}$, gravitational forces dominate vertical recharge to lower media without lateral discharge.

Figure 13d reveals decreasing discharge from PM I to the stream with reduced epikarst permeability. Cross-referencing Fig. 13a and e, lower epikarst permeability reduces both stream discharge and water level, limiting recharge to PM II. Figure 13f shows negligible epikarst permeability influence on PM III's hydrograph.

Epikarst permeability constitutes a critical factor in hydrological modeling of karst aquifer systems. Highly permeable epikarst produces rapid streamflow peaks followed by sharp declines, reflecting efficient groundwater leakage to the stream. Conversely, low permeability yields diminished peaks and broader discharge curves. While karst spring discharge remains relatively stable when epikarst permeability differs from porous media, proper characterization of epikarst permeability is essential for accurately simulating hydraulic interactions between media, regulating groundwater flow pathways and velocities. This enhances model reliability in capturing complex flow dynamics within karst conduit-stream systems.

4.3 Influence of Porosity on the Interaction Between Karst Aquifer Systems and Adjacent Streams

Figure 14 presents the hydrographic process curves under different porosity conditions (φ=0.4, φ=0.3, φ=0.2, φ=0.1) for (a) stream, (b) karst spring, (c) epikarst, (d) PM I, (e) PM II, and (f) PM III. Figure 14c.1 illustrates the schematic diagram of groundwater flow under different pore sizes. The study aims to elucidate how porosity regulates fluid flow patterns in complex conduit systems.

Figure 14 Hydrograph curves under different porosity conditions (φ=0.4, φ=0.3, φ=0.2, φ=0.1) for (a) stream, (b) karst spring, (c) epikarst, (d) PM I, (e) PM II, and (f) PM III. Among these, (c.1) illustrates a schematic diagram of the medium's water storage capacity and flow capacity under varying porosity conditions.

As shown in Fig. 14a, lower porosity results in higher flow peaks and earlier peak times. This occurs because reduced pore space limits groundwater storage capacity, forcing excess water to discharge rapidly and elevating the stream hydrograph. Figure 14b demonstrates that lower porosity drives groundwater to preferentially flow through karst conduits and discharge at springs. In Fig. 14c, the peak discharge of epikarst at φ=0.4 slightly exceeds those at φ=0.3, φ=0.2, and φ=0.1.

Figure 14d reveals that at φ=0.1, the storage capacity of PM I reach critical limits. Groundwater recharged from epikarst to PM I is rapidly discharged, resulting in significantly higher discharge rates compared to φ=0.3, φ=0.2, and φ=0.1. Figure 14e indicates increased discharge from porous media to the stream as porosity decreases. Combined with Fig. 14a, reduced porosity enhances stream stage and discharge but diminishes the stream's ability to recharge porous media due to limited storage capacity. Figure 14f shows negligible porosity effects on the hydrograph of PM III, as its behavior is primarily governed by conduit flow.

In hydrological modeling, porosity parameters must be calibrated to accurately simulate groundwater flow paths and storage-release dynamics. For low-porosity regions, models should emphasize rapid drainage capacity of conduit systems and transient flow variations. In high-porosity areas, considerations should include fluid retention risks, stream-porous media interactions, and their long-term impacts on geological stability and water resource allocation. Proper porosity parameterization enhances simulation accuracy for diverse hydrological processes, enabling improved prediction and management of karst water resources.

5 Conclusions

This study employed the Darcy–Brinkman–Stokes (DBS) method and the Volume of Fluid (VOF) technique to develop a unified model capable of coupling seepage and free flow, and meticulously characterizing two-phase (water-air) dynamics in a karst aquifer-stream system. The research confirms that, compared to conventional models like MODFLOW-CFPv2, this unified, multi-physics approach is essential for capturing the complex, dynamic processes inherent to karst systems.

High Non-linearity and Threshold Effects: The interaction between the karst aquifer and the stream is a highly non-linear process. Precipitation intensity acts as the primary driver, fundamentally altering flow paths and the contribution ratios of different media by triggering dynamic saturation, overflow, and synergistic recharge.

Necessity of Multi-Medium Coupling: The system's hydrological response is not governed by any single medium, but is co-determined by the rapid drainage capacity of conduits, the storage capacity of the matrix, and the permeability of the epikarst. For instance, while conduit geometry primarily controls peak discharge and recession efficiency, matrix porosity and epikarst permeability dictate the system's buffer capacity and the overall hydrograph morphology.

Importance of Unsaturated Zone Physics: The simulation results underscore the necessity of accurately describing unsaturated zone physics. The choice of Water Retention Models significantly impacts the stream hydrograph by altering the water storage and release dynamics of the matrix.

In summary, this study provides a robust framework for karst hydrological simulation. It demonstrates that a unified model capable of resolving coupled multi-medium and multi-phase flow is imperative for accurately predicting the complex hydrological responses of karst systems under varying precipitation scenarios. This enhanced predictive capability is fundamental for moving beyond oversimplified single-continuum models and developing more effective strategies for flood risk assessment, sustainable water resource allocation, and contamination vulnerability planning in these sensitive environments.

In future work, this research framework can provide critical tools for karst groundwater management:

By capturing non-linear thresholds, the model can more accurately predict how specific rainfall events trigger disproportionate flood peaks, thereby improving flood warning systems.

Aquifer Vulnerability Assessment: By coupling with a solute transport model, the framework can differentiate between acute/rapid contamination risks in conduits and chronic/slow risks in the matrix, providing a scientific basis for developing targeted source water protection strategies.