The article focuses on studying the interactions between hydrological processes involving karst aquifers, conduits, and streams, varying with precipitation and other factors that are later described in detail. The work is well-structured, the figures are of good quality, and the reading is fairly smooth.

 Dear Reviewer: Thank you very much for your comments and suggestions.  Our replies are listed as follow:

(1) The topic is, in my opinion, significant but somewhat self-contained when the discussion remains limited to numerical results without extending them to a broader analysis—namely, how this study could be useful in managing water resources in these complex aquifers or how it could be applied at a regional scale.

Thank you for pointing out the need to strengthen the practical application relevance in the Discussion section. In response to this suggestion, we have added implications of modeling parameters for karst water resources management in the revised version, as follows:

• Lines 847-855: Discusses the effect of conduit geometry (e.g., square cross-section) on aquifer-stream interaction fluxes, emphasizing its critical role in enhancing simulation accuracy.
• Lines 886-894: Analyzes the regulatory mechanism of epikarst permeability on groundwater flow paths, illustrating the importance of parameter calibration for water resources prediction.
• Lines 917-930: Adds an analysis of the impact of porosity parameters on karst water resources allocation and geological stability, clarifying the model's potential application in flood warning.

Lines 847-855:

“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.”

Lines 886-894:

“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.”

Lines 917-930:

“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.

Karst hydrological vulnerability manifests prominently through rapid infiltration, epikarst runoff, groundwater table fluctuations, and abrupt spring discharge variations. The DBS model effectively simulates multi-media interactions during extreme recharge events, enabling temporal analysis of media-stream exchanges, identification of peak interaction values, and applications in coupled conduit flow-seepage processes for two-phase flow systems.”

(2) The analytical and numerical approach is based on a representative sample of an aquifer, with the conduit size set by the authors (which I also find quite large) and a variation in rainfall intensity that I do not understand, as it is measured in meters according to Table 1. In hydrology, rainfall intensity is measured in mm/h, representing the amount of rainfall per unit of time. I kindly ask the authors to better define this aspect and, if necessary, correct it.

Thank you for your corrections regarding the units for conduit dimensions and precipitation intensity:

• Rationality of conduit dimensions: We have added a parameter sensitivity analysis (Lines 816-855), which quantifies the impact of different conduit diameters (0.2–0.5 m) and cross-sectional shapes (circular/square) on fluxes between media components using Fig. 12. This validates the physical basis for the parameter settings.
• Correction of precipitation intensity unit: The entry “Precipitation intensity (m)” in Table 1 was a typographical error; we have deleted this column. Furthermore, detailed definitions of the precipitation function have been added in the “5 Rainfall Infiltration Recharge Boundary” subsection (Lines 712-735): The precipitation function is defined as a time-dependent variable I(t), with its intensity modulated by the dimensionless parameter b, consistent with conventional expressions in the hydrology field.

Lines 816-855:

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

Fig. 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) porous medium I (PM I), (e) PM II, and (f) PM III. Fig. 12(c.1) illustrates different conduit cross-sectional shapes to analyze their impacts on the interactive flow between karst aquifer systems and adjacent streams.

As shown in Fig. 12(a), 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_c=0.5) exhibits higher peak discharge than its circular counterpart (rc=0.5) due to its surplus cross-sectional area accommodating greater fluid discharge under identical nominal radii.

Fig. 12(b) demonstrates that karst spring peak discharge increases with conduit radius. At r=0.5 m, the square-section conduit (s-r_c=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 Figs. 12(c), (d), and (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.

Fig. 12(e) indicates that larger conduit radii correspond to lower negative values. Combined with Fig. 12(a), 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. 12(f) 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.”

Lines 712-735:

“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):

 Here,  denotes the time of the -th rainfall event, and  represents the total rainfall amount at that time. According to Chang et al. (2015), the parameters 、 , and  are set as constants (90, 1.5, and 20, respectively). Variations in rainfall intensity during the infiltration recharge process are controlled by adjusting the value of the dimensionless parameter .”

(3) Another aspect concerns the discussion section. This should compare the obtained results with those in the literature, discuss the limitations and advantages of such an approach, perform a sensitivity analysis on the calibrated parameters, and assess the study’s usefulness and general applicability. The current section focuses on comparing the results with those from the MODFLOW-CFP model, which the authors themselves implemented. In my opinion, this should still be described in the results section, and consequently, the methodological section should be expanded to include this additional approach.

Thank you for your suggestions regarding results comparison and chapter organization:

• Literature comparison: Added quantitative validation against experimental data from Warrick et al. (1985) and Vauclin et al. (1979) (Lines 712-735), demonstrating the reliability of the DBS model in simulating variably saturated flow.

• Chapter restructuring: Relocated the MODFLOW-CFPv2 comparison section to the "Results" chapter (Lines 352-354) and expanded the explanation of CFPv2 principles in the "Methodology" section (Lines 328-340).

Lines 712-735:

“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.

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:  = 1 × 10−6 m/s,  = 0.363,  = 0.186, and n = 1.53. The model consists of a vertical soil column (1  thickness) with an initial pressure head of -8  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  in height, 6.00  in width, and 0.05  in thickness, with an impermeable base and free drainage boundaries on both sides. Initially, the water table is set at 0.65 . A central 1.00  section of the top boundary receives uniform precipitation at 0.148 m/h for 8 hours, during which free surface evolution is monitored. Soil hydraulic properties are described using the van Genuchten-Mualem model with parameters:  = 0.35 ,   = 0.30,  = 0.01. Due to symmetry, the DBS model simulates the right half of the domain.”

Lines 328-340:

“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. (2020), Gao et al. (2020), and Shirafkan et al. (2023).”

(4) line 165: wrong unit of measure for permeability, gravitational acceleration

Thank you for your comment. We have corrected the labeling errors in the revised manuscript.

(5) line 205: delete repetitive title.

Thank you for your comment. We have removed the duplicate headings in the revised manuscript.

(6) Figure 1: insert letters a and b also in the figure for better readability.

Thank you for the suggestion. We have clearly labeled subplots (a), (b), (a.1), and (b.1) in Fig. 1 and updated the captions (Lines 1200–1207) to enhance readability.

General Comments:

This manuscript presents a coupled Darcy-Brinkman-Stokes (DBS) model to simulate karst aquifer and stream interactions under various precipitation conditions, comparing its results to MODFLOW-CFP. The study is technically sound, and the modeling framework has clear potential for advancing process-based simulation of turbulent flows and unsaturated conditions in karst systems. However, the manuscript currently suffers from significant shortcomings in clarity, structure, model transparency, and broader contextual discussion. These issues limit its accessibility and impact.

Thank you for your recognition of the model framework and constructive feedback on improving the manuscript. We have thoroughly revised the content, with specific details outlined as follows:

(1)  A clearer articulation of the research gap and contributions in the Introduction.

Sincere thanks for your comments. We have supplemented the research gap and this study’s contributions in the Introduction section:

·       Clarified research gaps: Noted limitations of CFPv2 (single-phase flow only) and Hydrus (difficulty in coupling conduit flow) (Lines 92–94, 106–108), highlighting the innovation of the DBS model in fully coupling two-phase flow, variably saturated flow, and turbulent flow (Lines 116–123).

·       Enhanced contribution statements: Added a concluding paragraph (Lines 132–139) emphasizing the model’s breakthrough in resolving karst aquifer-stream combined recharge mechanisms.

(2)  A reorganization and simplification of the Methods section, with better-defined subsections and clearly stated variable definitions.

Thank you for your comments. We have revised the Methodology section and provided a variable definition table.

• Structural adjustments: Reorganized into subsections 2.1–2.5 (including the variable definition table, lines 1191-1193) to enhance logical clarity.

• Clarification of key concepts: Added justification for turbulence simulation (lines 256-260), definition of two-phase flow parameters (lines 211-237), and the rationale for selecting the VGM/BC models (lines 294-326).

(3) Improved figure quality and more structured result narratives that highlight key trends, comparisons, and implications.

Thank you for your comments. In the revised manuscript, we have reorganized the results section and included a comparison of CFPv2 with our study.

(4) A more rigorous and quantitative comparison with MODFLOW-CFP, beyond qualitative differences.

Thank you for your comments. In the revised manuscript, we have conducted a quantitative comparison between CFPv2 and the DBS study, and provided a quantitative comparison table.

•   Quantitative comparison with CFPv2: Added Table 3 (line 1197) comparing metrics such as peak discharge and peak arrival time to clarify the advantages of the DBS in capturing double flood peaks.

(5) A deeper discussion of how the findings contribute to our understanding of karst-stream interactions and why this modeling advance matters in practice.

Specific Comments Title

(6) The current title is vague and lacks specificity about the method or novelty.

• Consider revising to reflect the modeling framework or scientific contribution, e.g.,

“Simulating Precipitation-Induced Karst-Stream Interactions Using a Coupled Darcy-Brinkman-Stokes Model”.

Thank you for your comments. We have revised the title in the updated manuscript.

• Lines 1-2:“Simulating Precipitation-Induced Karst-Stream Interactions Using a Coupled Darcy-Brinkman-Stokes Model ”

Abstract:

(7) Long abstract detracts from its purpose, it should be short and to the point.

In the revised manuscript, we have condensed the abstract length and clarified the main focus of the paper.

Lines 14-45:

“Abstract: The variation in seasonal precipitation intensity impacts the dynamic interaction between the karst aquifer and stream.  However, the interaction mechanism between the karst aquifer and stream is currently still unclear, and characterizing the impact of dynamic saturation process of groundwater in karst media on the interaction process remains a challenge.  This study provides an in-depth analysis of the interaction processes between karst aquifer systems and adjacent streams, along with water-air two-phase flow in aquifer media. Multiple water retention models were employed to characterize the soil-water characteristics of porous media and variably saturated groundwater flow. The research reveals that rainfall intensity variations significantly influence the interactions between karst aquifer systems and streams. These interactive processes become increasingly complex with higher rainfall intensities, involving multi-media collaborative recharge and dynamic interactions, while the contribution proportions of different media to streamflow also change accordingly. By comparing the modeling differences and numerical results between CFPv2 and DBS approaches in generalized models, the validity of the DBS model for groundwater modeling was verified. Under consecutive rainfall events, total rainfall intensity plays a crucial role in hydrological process variations of adjacent streams. Groundwater stored in porous media of karst systems during the first rainfall event was found to influence stream water levels during subsequent rainfall events, while conduit storage exhibited minimal impact. Multi-level conduit configurations under specific conditions, particularly during intense rainfall, can significantly affect hydrological processes in both streams and karst conduits. Uncertainty analysis demonstrates that conduit geometry, diameter, epikarst permeability, and porosity differentially influence hydrological processes in karst aquifer systems. Variations in these parameters induce corresponding changes in peak flow rates, peak timing of stream and karst spring discharges, as well as redistribution of discharge contributions among different media, ultimately affecting the overall hydrological dynamics of the coupled karst aquifer-stream system.  It can accurately depict the two-phase interactive flow between various media controlled by the dynamic saturation process, and reveal the dynamic interaction process between karst aquifers affected by the epikarst, sinkholes, and conduits under infiltration recharge and stream.  Meanwhile, it can precisely explain the processes of infiltration, overflow, and recession.”

(8) Avoid referring to equation names (e.g., DBS, BC, VGM) without explanation. Focus on what was learned, not just what was done.

Thank you for your comments. We have rewritten the abstract in the revised manuscript to avoid direct reference to equations.

Introduction:

(9) Introduce a clear knowledge gap: why are existing models insufficient?

In the revised manuscript, we have clarified the rationale for employing the Darcy-Brinkman-Stokes (DBS) model based on the following key points:

• Applicability Limitations of MODFLOW-CFP: The applicability of MODFLOW-CFP is limited by its current capability to simulate only single-phase flow (Lines 92-94).
• Inherent Difficulties with Hydrus: The Hydrus model faces inherent difficulties in coupling rapid conduit free flow with porous media seepage flow in karst regions, due to its lack of a specialized conduit flow solution (Lines 106-108).
• Advantages of the DBS Model: In contrast, the DBS model emerges as an ideal choice for simulating the complex flow within karst systems (including saturated-unsaturated two-phase flow) because of its theoretical foundation for coupling free flow and seepage flow within a unified framework. Multiple recent studies cited in Lines 116-123 (Huang et al., 2024; Nillama et al., 2022; Carrillo et al., 2020; Lu et al., 2023; Soulaine, 2024) adequately demonstrate the effectiveness and potential of the DBS approach in characterizing the coupled processes of conduit/fracture flow and matrix seepage flow.

Lines 92-94:

“Although MODFLOW-CFP is relatively comprehensive for regional karst groundwater simulation studies, the current version of MODFLOW-CFP only supports modeling single-phase groundwater flow.”

Lines 106-108:

“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.”

Lines 116-123:

“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.”

(10) State specific objectives and contributions at the end of the introduction.

We added the following statement at the end of the Introduction (Lines 132-139):

The DBS model represents the first implementation enabling the fully coupled simulation of two-phase flow, variably saturated conditions, and turbulent flow within karst systems. This addresses a gap in existing tools (e.g., CFPv2), which cannot characterize the dynamic, synchronized recharge processes involving saturation changes.

Lines 132-139:

“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.”

(11) Consider breaking the introduction into shorter paragraphs or subsections.

Thank you for the comment. In the revised manuscript, we have divided the longer paragraphs to enhance readability.

Methods:

(12) The modeling section is too dense and difficult to follow.

We appreciate your valuable feedback. In the revised manuscript, we have thoroughly restructured the Methods section as follows:

2.Materials and methods

2.1 Numerical modelling

2.2 DBS model

2.2.1 Two-Phase Flow Parameter Definition

2.2.2 Governing Equations

2.2.3 Subdomain Formulation

2.2.4 Relative Permeability Model

2.3 CFPv2 model

• Model Comparison and Numerical Model Construction

2.4.1 DBS Model Conversion and Applicability Assessment

2.4.2 Model Comparison and Discretization Schemes

2.5 Rainfall Infiltration Recharge Boundary

(13) Key concepts (e.g., turbulent conduit flow, two-phase flow, saturation) are introduced without sufficient explanation.

Thank you for the comment. In the revised manuscript, we have provided clearer explanations of turbulent conduit flow, two-phase flow, and saturation dynamics before reintroducing these concepts.

Lines 256-273:

“Conduit networks in karst aquifer systems are often associated with turbulent flow (Reimann et al., 2011). To resolve turbulence in the DBS (Dual-domain Brinkman-Stokes) equations, the Unsteady Reynolds-Averaged Navier-Stokes (URANS) framework is required. As demonstrated by del Jesus et al. (2012), the k-epsilon turbulence model is effective for evaluating turbulent processes within porous media. Consequently, the k-epsilon-based DBS turbulence governing equations are formulated as follows:

where,  represents the turbulent velocity vector [ ],  is the relative velocity of gas-phase and water-phase turbulence [ ], and  is the effective viscosity, which can be defined as , where  is the dynamic viscosity and  is the turbulent kinetic viscosity.

The eddy viscosity is expressed as:           

where: ​: Turbulent kinetic energy per unit mass [ ], : Turbulent dissipation rate per unit mass [ ], ​: Dimensionless constant with a value of 0.09.”

Lines 211-237:

“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.

In this context,  represents porosity,  denotes the total volume of the unit [ , while and ​ correspond to the volumes of the liquid phase (water) and gas phase (air), respectively [ .

Hirt and Nichols (1981) introduced the Volume of Fluid (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 , where:

Liquid phase saturation: .

Gas phase saturation: .

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 :

The average fluid density  and viscosity  [ ] within a grid cell are calculated via saturation-weighted averaging:

where ​ is the gas phase density  and  is the liquid phase (water) density .

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

where: is the fluid velocity vector [ ],  is the relative velocity between the gas and liquid phases [ ].“

(14) Variables are often undefined or defined long after their use.

We appreciate your comment. In the revised manuscript, we have carefully refined the definitions of key variables to improve clarity and precision.

(15) Provide a table of variable definitions.

We appreciate your comment. In the revised manuscript, we have included a comprehensive variable definition table (see Table 1) to clarify all key parameters and their respective units.

(16) Clearly justify modeling choices: Why use turbulence? Why VGM vs. BC?

Thank you for the comment. In the revised manuscript:

• Regarding turbulence justification
  
  We have clarified the rationale for selecting turbulent flow modeling (Lines 256-260).
• Concerning VGM vs. BC model comparison
  
  We now provide a parallel analysis of both the VGM and BC models, with specific discussion of scenarios where the VGM formulation demonstrates advantages over the BC approach (Lines 294-326).

Lines 256-260:

“Conduit networks in karst aquifer systems are often associated with turbulent flow (Reimann et al., 2011). To resolve turbulence in the DBS (Dual-domain Brinkman-Stokes) equations, the Unsteady Reynolds-Averaged Navier-Stokes (URANS) framework is required. As demonstrated by del Jesus et al. (2012), the k-epsilon turbulence model is effective for evaluating turbulent processes within porous media.“

Lines 294-326:

“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:

where: denotes the effective water saturation,  and  represent the water saturation and residual water saturation, respectively, and  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:

denotes the relative permeability, where  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 & Or, 2013). The van Genuchten (1980) model addresses this limitation. By incorporating the parameter  proposed by Mualem (1976), the modified van Genuchten-Mualem (VGM) model (Parker et al., 1987) is formulated as:

Here, m is a dimensionless parameter.

The selection of permeability equations is critical for appropriate predictions of relative permeability (Yang et al., 2019), 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.”

Results:

(17)The results are structured around simulations but lack comparative synthesis. Each simulation is described in isolation.

Thank you for the comment. In the revised manuscript, we have restructured the Results section according to the following framework:

1. Results

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

3.1.1 Karst Aquifer-Stream Interactions Under Varying Precipitation Intensities

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

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

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

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

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

(18) Start each subsection with a motivating question (e.g., "How does conduit depth affect stream discharge?").

We appreciate your insightful observation regarding the value of engaging questions. In the revised manuscript, we have implemented the following enhancements to improve reader engagement:

Lines 557-558:

“How does the threshold attainment of storage capacity in the lower porous media affect the hydrological processes of the upper porous media ?”

Lines 633-636:

“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 ?”

(19) Emphasize key differences in system behavior (e.g., thresholds, flow reversals, saturation delays).

Thank you for your comments. We have actively incorporated your suggestions in the revised manuscript to analyze key systematic differences.

(20) Integrate comparisons to MODFLOW-CFP throughout instead of relegating to the end.

We have conducted a comparative analysis between DBS and CFPv2 in the section preceding the two-phase flow results.

Comparison with MODFLOW-CFP:

Lines 478-494:

“Based on the comparison between DBS and Modflow-CFPv2 results in Figs 4(a), (b), and (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.“

(21)The comparison is largely qualitative.

Thank you for your comments. We have added quantitative comparisons in the revised manuscript.

• Add quantitative metrics (e.g., time to saturation, max stream recharge, total outflow).

We have supplemented relevant quantitative metrics in the revised manuscript. Peak Lag Time ；Peak Flow ；Total Outflow

• Include a summary table comparing DBS and MODFLOW results for key variables.

Added Table 3 (Line 1197) with quantitative comparisons of key metrics (time-to-peak, peak discharge, total outflow), demonstrating that the DBS model effectively captures complex processes such as secondary flood peaks (which CFPv2 overlooks due to its simplified vertical layering). This highlights DBS's superior capability in characterizing multi-media interactions.

(22) Discuss the limitations of MODFLOW-CFP more directly — why does it fail to simulate certain behaviors ?

Thank you for your comments.

• The DBS model successfully captured the secondary flood peak overlooked by CFPv2, revealing a multi-media synergistic recharge mechanism.
• CFPv2 failed to simulate the overflow-driven secondary flood peak (Fig. 3) due to its simplified vertical layering (1D conduit flow), whereas DBS overcame this limitation through coupled 3D conduit flow modeling.

Lines 352-354:

“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.“

Discussion:

(23) The discussion largely restates results and does not engage deeply with broader hydrologic or modeling implications.

Thank you for your comments. We have reorganized the results in the Discussion section and further investigated the influence of sensitivity parameters on the modeling:

• Address: What do these findings imply about karst vulnerability under extreme precipitation?
• Discuss potential applications of the DBS model in field calibration, water management, or integrated modeling platforms.

Thank you for pointing out the need to strengthen the connection to practical applications in the Discussion section. In the revised manuscript, we have supplemented the implications of modeling parameters for water resource management in karst areas, including the vulnerability of karst systems under extreme precipitation events using the DBS model, and its potential applications in water resource management or integrated modeling platforms:

• Lines 847–855: Discusses the impact of conduit geometry (e.g., square cross-section) on aquifer-stream exchange fluxes, emphasizing its key role in enhancing simulation accuracy.
• Lines 886–894: Analyzes how epikarst permeability governs groundwater flow paths, illustrating the importance of parameter calibration for water resource predictions.
• Lines 917–930: Adds an analysis of the influence of porosity parameters on karst water resource allocation and geological stability, clarifying the model's potential applications in flood early warning systems.

Lines 847–855:

“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.”

Lines 886–894:

“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.”

Lines 917–930:

“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.

Karst hydrological vulnerability manifests prominently through rapid infiltration, epikarst runoff, groundwater table fluctuations, and abrupt spring discharge variations. The DBS model effectively simulates multi-media interactions during extreme recharge events, enabling temporal analysis of media-stream exchanges, identification of peak interaction values, and applications in coupled conduit flow-seepage processes for two-phase flow systems.”

(24) Compare to other recent modeling studies in karst hydrology (e.g., those using HYDRUS, ParFlow, or CFP enhancements).

Thank you for your comments. We first replaced the comparative CFP model with the updated CFPv2 model and added two references for comparative validation.

• Literature Comparison: Added quantitative validation against experimental data from Warrick et al. (1985) and Vauclin et al. (1979) (lines 712-735), demonstrating the reliability of the DBS model in simulating variably saturated flow.
• Model Comparison: Compared the DBS model with the enhanced CFP version (MODFLOW-CFPv2) (lines 352-354), and expanded the explanation of CFPv2 principles in the "Methods" section (lines 328-340). (29)

Missing Elements:

(25) Model Validation: No benchmark or field validation is presented — even conceptual model validation would strengthen the work.

Thank you for your comments. We conducted comparative validation of the model using two published studies.

• Literature Comparison: Added quantitative validation against experimental data from Warrick et al. (1985) and Vauclin et al. (1979) (lines 712-735), demonstrating the reliability of the DBS model in simulating variably saturated flow.

(26) Uncertainty Analysis: Sensitivity to soil properties, precipitation rate, conduit geometry, etc., is not discussed.

Thank you for your comments. The revised manuscript discusses sensitivity analyses and potential impacts of soil properties, conduit size, conduit shape, porosity, and permeability:

• A new “4. Uncertainty Analysis” section (lines 810-930) demonstrates that:

(1) Conduit geometry (square sections yield higher peak flows than circular ones)

(2) Permeability of the epikarst zone (higher permeability enhances inter-media exchange frequency)

are key factors regulating karst-stream feedback mechanisms. This provides prioritization guidance for model parameter calibration.

1. 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

Fig. 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) porous medium I (PM I), (e) PM II, and (f) PM III. Fig. 12(c.1) illustrates different conduit cross-sectional shapes to analyze their impacts on the interactive flow between karst aquifer systems and adjacent streams.

As shown in Fig. 12(a), 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_c=0.5) exhibits higher peak discharge than its circular counterpart (rc=0.5) due to its surplus cross-sectional area accommodating greater fluid discharge under identical nominal radii.

Fig. 12(b) demonstrates that karst spring peak discharge increases with conduit radius. At r=0.5 m, the square-section conduit (s-r_c=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 Figs. 12(c), (d), and (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.

Fig. 12(e) indicates that larger conduit radii correspond to lower negative values. Combined with Fig. 12(a), 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. 12(f) 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. Fig. 13 illustrates the hydrological process curves under different epikarst permeability coefficients (K_E=10⁻⁶, 10⁻⁷, 10⁻⁸, 10⁻⁹; when K_E=10⁻⁹, 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.

As shown in Fig. 13(a), under high epikarst permeability (K_E=10⁻⁶): the discharge curve rises rapidly to a peak of ~4.5  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.

Fig. 13(b) demonstrates that epikarst permeability differences from porous media have minimal impact on conduit flow. However, when epikarst permeability equals that of porous media (K_E=10⁻⁹), the peak discharge at the karst spring decreases while maintaining identical baseflow recession characteristics. Combining Figs. 13(c) and (c.1), higher epikarst permeability enhances lateral discharge to the stream. At K_E=10⁻⁹, gravitational forces dominate vertical recharge to lower media without lateral discharge.

Fig. 13(d) reveals decreasing discharge from Porous Medium I to the stream with reduced epikarst permeability. Cross-referencing Figs. 13(a) and (e), lower epikarst permeability reduces both stream discharge and water level, limiting recharge to Porous Medium II. Fig. 13(f) shows negligible epikarst permeability influence on Porous Medium 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

Fig. 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. Fig. 14(c.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.

As shown in Fig. 14(a), 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. Fig. 14(b) demonstrates that lower porosity drives groundwater to preferentially flow through karst conduits and discharge at springs. In Fig. 14(c), the peak discharge of epikarst at  =0.4 slightly exceeds those at  =0.3,  =0.2, and  =0.1.

Fig. 14(d) reveals that at  =0.1, the storage capacity of porous medium I reaches critical limits. Groundwater recharged from epikarst to porous medium I is rapidly discharged, resulting in significantly higher discharge rates compared to  =0.3,  =0.2, and  =0.1. Fig. 14(e) indicates increased discharge from porous media to the stream as porosity decreases. Combined with Fig. 14(a), reduced porosity enhances stream stage and discharge but diminishes the stream's ability to recharge porous media due to limited storage capacity. Fig. 14(f) shows negligible porosity effects on the hydrograph of porous medium 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.

Karst hydrological vulnerability manifests prominently through rapid infiltration, epikarst runoff, groundwater table fluctuations, and abrupt spring discharge variations. The DBS model effectively simulates multi-media interactions during extreme recharge events, enabling temporal analysis of media-stream exchanges, identification of peak interaction values, and applications in coupled conduit flow-seepage processes for two-phase flow systems.”