In karst aquifers, groundwater flow is heavily influenced by the interconnected underground cavities and conduits known as karst networks. Spatial modeling of karst flows requires the use of distributed approaches accounting for these networks. Due to the highly variable conduit thickness, their exploration is, however, often complex, and mapping them using indirect methods such geophysical methods as seismic reflection or electrical resistivity has proven challenging. To address these limitations, stochastically simulating discrete karst networks permits to account for the uncertainties on conduit position and geometry. Among the existing methods proposed in the literature, only a few can reproduce realistic and diverse karst morphologies. In this work, we propose a new approach for simulating discrete karst networks, that incorporates field data to generate a range of possible models. It solves a shortest path problem between sinks and springs – respectively the inlets and outlets of the network – with the use of an anisotropic cost function defined on an unstructured mesh conformal to geological and structural heterogeneities. This cost function represents the physico-chemical processes that govern speleogenesis – such as erosion and chemical weathering – providing simplified control over the morphometry of the generated networks. Our approach reproduces the vadose-phreatic partition visible in the karst networks, by generating subvertical conduits in the unsaturated zone and subhorizontal ones in the saturated part. It encompasses geological parameters such as inception surfaces, fractures, permeability, and solubility of layers, along with considering the hydrological context of recharge by assigning relative weights to the inlets. To evaluate our approach, we simulate various synthetic models to demonstrate the influence of different input parameters on the spatial organization of past and present karst flows. We also apply our method to a real case study : the Fontaine de Ribeaucourt (Meuse, France).