In karst aquifers, groundwater flow is highly influenced by the interconnected underground cavities and conduits that form the karst network. Modeling karst flows requires the use of spatially distributed approaches accounting for these networks. Their exploration is, however, often complex, and mapping them using indirect methods such as geophysical ones has proven challenging. To overcome these limitations, stochastically simulating discrete karst networks should account for the uncertainties on conduit position and geometry. Only a few existing methods can reproduce realistic and diverse karst morphologies. In this work, we propose a new algorithm, KarstNSim, for simulating discrete karst networks, that incorporates field data to generate a range of possible karst network geometries. It relies on the computation of the shortest path between the inlets and outlets of the network with the use of an anisotropic cost function defined on an n-nearest neighbor graph 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 geometry of the generated networks. Our approach reproduces the vadose-phreatic partition visible in the karst networks, by generating sub-vertical conduits in the unsaturated zone and sub-horizontal 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. We simulate various synthetic models to demonstrate the influence of different input parameters on the spatial organization of the past and present karst flows. We also apply KarstNSim to a real case study-the Ribeaucourt system (Meuse, France)-and use available data to generate realistic karst networks. The lack of knowledge on the Ribeaucourt network geometry is alleviated by the use of multiple scenarios corresponding to different sets of weights in the cost function, resulting in a variety of network morphologies which represent the uncertainty on the real geometry. The geometrical and topological metrics computed on the simulated networks generally overlap with those of real ones, suggesting the reliability of the method.