Numerical Simulation Study of Preferential Flow at the Rock-Soil Interface in the Epikarst Zone

In the karst region of southwest China, the epikarst zone is characterized by soil-filled solution grooves and an intricate network of fractures. Rainfall infiltration is jointly controlled by the groove soil, preferential flow pathways at the rock-soil interface, and the surrounding fractures. However, the dynamic interactions among the groove soil, rock-soil interface, and fracture network, especially their response mechanisms under different rainfall intensities, remain poorly quantified. Based on in-situ dye tracer tests on a solution groove profile, combined with long-term monitoring of soil moisture and fracture water content, this study developed a numerical model of water flow to quantitatively reveal the contributions of preferential flow at the rock-soil interface and adjacent inclined fractures to water movement and the outlet flux at the lower boundary. The study was conducted on a typical epikarst profile in Puding, Guizhou Province. The profile is 10 m long and 2.24-2.5 m high, containing three solution grooves with widths of 0.22-0.82 m, four sub-horizontal fractures (dip angle~5.4°) and three sub-vertical fractures (dip angle~83.2°). A dye tracer test using Brilliant Blue FCF was carried out on the middle solution groove at a rainfall intensity of 60 mm·h⁻¹ for one hour, followed by excavation to reveal the spatial distribution of preferential flow. Five soil moisture sensors were installed at depths of 5, 20, 50, 100 and 180 cm in the groove soil, and five fracture moisture sensors (F 1-F 5) were placed in adjacent fractures. Meteorological data, including rainfall, were recorded every 30 minutes from September 2022 to September 2024. Based on the soil and fracture moisture responses, four representative rainfall events were selected for detailed analysis: light rain (4.3 mm), moderate rain (12.3 mm), heavy rain (34.2 mm), and storm (62.1 mm).A numerical model of water flow was constructed using COMSOL Multiphysics. Unsaturated flow in the groove soil and fractures was described by the Richards equation, and the soil water retention and hydraulic conductivity were described by the van Genuchten model. Solute transport in the dye tracer experiment was simulated using the convection-diffusion equation with Langmuir isothermal adsorption. The effective aperture of the open rock-soil interface was inversely calibrated against the dye distribution and soil moisture dynamics, yielding an aperture of 4.22 mm in the upper zone (A) and 1.66 mm in the lower zone (B). Hydraulic parameters for the soil and fractures were optimized using the BOBYQA iterative algorithm. Model validation against two years of soil and fracture moisture data gave root-mean-square errors (RMSE) mostly below 0.05 m³·m⁻³ and mean absolute percentage errors (MAPE) below 5%, indicating that the calibrated parameters represent the average hydraulic conditions over the entire monitoring period.Three simulation scenarios were compared under a storm event: (A) actual conditions with an open rock-soil interface (preferential flow) and surrounding fractures; (B) closed rock-soil interface (no preferential flow); and (C) no surrounding fractures. The results show that preferential flow at the rock-soil interface provides an efficient pathway for vertical water movement. Under actual conditions (scenario A), the maximum flow velocity reaches 0.027 m·s⁻¹. When preferential flow is ignored (scenario B), the flow exhibits typical piston-flow characteristics with a uniform flow field and a maximum velocity of only 7.32$ \times $10⁻⁵ m·s⁻¹. When surrounding fractures are absent (scenario C), the flow velocity in the groove is about 1.3 times that of scenario A, indicating that inclined fractures retard vertical flow under heavy rainfall. For the four representative rainfall events, the outlet water flux at the lower boundary of the solution groove was analyzed. In the preferential flow zone of the rock-soil interface, the flux increased significantly with rainfall intensity. In the non-preferential flow zone (soil matrix), the flux was very small or even negative during light and moderate rain events, indicating that infiltrated water did not reach the lower boundary; only after heavy rain and storm events did a slow increase occur, reflecting slow matrix flow. The cumulative outflow volumes over 24 hours were 0.19 L (light rain), 18.25 L (moderate rain), 29.62 L (heavy rain), and 132.33 L (storm rain) in the preferential flow zone, whereas the non-preferential zone contributed only 0.12-0.23 L for all events. This demonstrates that the preferential flow zone is one to two orders of magnitude more effective in recharging the deep epikarst than the non-preferential zone. Furthermore, the presence of inclined fractures significantly influences the outflow process. In the preferential flow zone, after heavy rain and storm events, ignoring the surrounding fractures led to a faster rise and a larger flux, indicating that the inclined fractures weaken the preferential outflow. In the non-preferential flow zone, after heavy rain and storm events, the increase in outlet flux when neglecting the surrounding fractures was generally greater than when fractures were present, suggesting that fractures also modulate matrix flow, but with a more persistent effect. In conclusion, this study confirms that preferential flow at the rock-soil interface is the dominant pathway for rapid groundwater recharge from precipitation. The combination of dye tracer tests, long-term in-situ monitoring, and physically-based numerical modeling provides a reliable approach for quantifying the roles of rock-soil interface aperture and fractures under different rainfall intensities. The findings offer a scientific basis for the coordinated management of water resources and ecological protection in karst regions.