The characterization and management of karst systems is a challenging task due to their inherently heterogeneous nature and vulnerability with respect to contamination. Highly conductive features of the vadose zone (e.g., dissolution shafts and faults) induce flow channeling and preferential flow. This complicates any efforts to simulate rapid recharge dynamics in deep porous-fractured vadose zones in the context of flood and contamination risk assessment. Therefore, a strong need for numerical modeling strategies arises that employ conceptually sound formulations of these dynamics based on physical processes. Here, we present a novel modeling strategy by extending the numerical discrete-continuum flow model MODFLOW-CFPv2 to allow the simultaneous computation of diffuse fluxes and film-flow in the vadose zone, thus simulating infiltration via preferential pathways. We conduct a global sensitivity analysis of a synthetic karst system that addresses the importance of including such processes in karst modeling. While event-averaged sensitivities are in alignment with commonly observed dominance of the phreatic zone properties, results of time-dependent sensitivities suggest that during strong infiltration events the consideration of film-flow and its controlling parameters, that is, the fracture facial area density and an applied upper threshold for its activation, can become important. Our distributed numerical method assists in the development of karst modeling strategies where a sufficiently large and developed vadose zone offers the capacity for preferential flow that may not be accurately reproduced by most bulk-effective methods. Hence, it benefits the unique characterization of such systems and can be easily implemented in existing workflows such as CFPv2.