This work presents a coupled numerical model to study the surface propagation of flows generated by geysering effects in various scenarios, including urban areas and erodible surfaces. The model combines a numerical-analytical approach for the formation and propagation of axisymmetric cavities with a reduced-depth model incorporating air-water interfaces. This coupling enables a comprehensive analysis of both the formation and propagation in conduits and the surface flow evolution and erosion over initially dry surfaces, considering both non-erodible and erodible media. A particular feature of this work is the use of the non-oscillatory finite element method NFEM. This numerical approach is a conservative sign-preserving algorithm that ensures an accurate representation of evolving interfaces and wavefronts, minimizing spurious oscillations and keeping physical coherence in the air-water-sediment transitions. This is particularly relevant for predicting complex phenomena associated with mixed flows in dynamically changing terrains. The proposed model represents a robust tool for simulating and predicting flows driven by geysering. It enables the study of propagation dynamics not only in urban areas but also in natural terrains susceptible to erosion.