This paper describes a methodology to generate stochastic models of 3D cave systems. It combines existing geostatistical techniques, namely object-based and variogram-based geostatistics. The first is used to simulate Discrete Fracture Networks (DFNs) which in addition of bedding planes drive the development of karst features. The discontinuity network and the matrix are jointly discretized into a graph of connectivities. Preferential flow paths determining the topology of the cave system are then extracted using graph search algorithms. The distance map to the selected paths is used to conditioned the spatial extent of karst features, i.e., the probability to simulate karst decreases moving away from “dissolution paths”. Eventually, multiple realizations of a distance cutoff are generated with Sequential Gaussian Simulation to perturb the distance function. Resulting models efficiently reproduce the spatial variability of karst features in term of geometry and their spatial organization in term of connectivity.