The numerical model proposed in this paper represents a new method in predicting initial karst spatial distribution of carbonate aquifers. Karstification Index, \( \varPhi_{\text{K}} \), is a measure of the possible largest amount of calcite that can be dissolved from a unit area of rock and transported by groundwater flow during a unit length of time. Mathematically, it is defined as \( \varPhi_{\text{K}} = q \cdot s \), where q is the flux of groundwater and s is the rock solubility; and Karstification Index, \( \varPhi_{\text{K}} \), is the function of flow, heat, and calcium concentration. The numerical model simulates the time- and space-dependent susceptibility to, and process of, karstification by iterative calculations between \( \varPhi_{\text{K}} \) and hydraulic conductivity. A sample calculation is presented on a hypothetical limestone plateau with steep escarpment, a spring, and a river on one side, a vertical no flow boundary on the other side, and a horizontal impermeable base. The rock has an equivalent porous-medium hydraulic conductivity of K = 5E−8 ms−1, which is characteristic to non-karstified limestone. Initially, the spring is located at the intersection of the water table and the escarpment’s face some distance above the river. The simulation result showed that the water table is lowered in 2325 years across the plateau, but the spring’s drop is stopped after 65 years. Karstification is most intensive along the water table due to the higher value of Karstification Index which resulted from the intensive flux of the water. The continual decline of the water table results in a laterally elongated and progressively expanding body of karstified rock and finally as a horizontal conduit at river level. The concentration of Ca2+ and conduit diameter decrease from the spring toward the impermeable vertical wall on the opposite side and downward to the base of the plateau. The results are consistent with general theories and empirical observations of karst and cave evolution.