We develop a new numerical model that couples hydraulic, mechanical and chemical processes to investigate the dynamic influence of stress on incipient karst formation in natural fracture networks. Our geological models are constructed based on realistic fracture patterns mapped from outcrops, which capture a wide range of spatial distribution and organization of carbonate fracture networks. We simulate the dissolutional growth of karst conduits in the fractured formation under varied initial aperture distribution and in-situ stress loading. We found that the importance of the stress effect on karstification depends on the relative relationship between the flow direction and structural anisotropy of the fracture network. When the flow occurs in the direction of main network structures, karst conduits only develop locally along a few large fractures with a preferential orientation for frictional sliding under the differential stress due to enhanced transmissivity caused by the important shear-induced dilation. In contrast, when flow is in the direction transverse to the main fractures, the far-field stress loading has a negligible impact on the emergent dissolution pattern while only impact on the onset time of breakthrough. In this case, the developed conduits are much more tortuous with many branches. In both cases, the presence of initial aperture variability facilitates the development of a dissolution front and therefore makes breakthrough faster. Our results demonstrate that the heterogeneity induced by geometrical complexities and in-situ stress conditions may play a decisive role in the karstification processes in fractured rocks. The results from this research provide important insights into the spatial relationship between tectonic structures and karst cavities, which has important implications for karst aquifer management and related environmental protection, or hydrocarbon production in fractured carbonate reservoirs.