<p>Zonal disintegration is a typical bifurcation phenomenon observed in the deep rock mass. A hardening-softening model was developed to investigate this failure behavior. In addition, the strength heterogeneity within the rock mass was incorporated via three discrete random distributions. The proposed approach was implemented and coupled with the finite difference method. Zonal disintegration was successfully captured in the numerical simulations. A strong correlation was identified between weak zones and the occurrence of zonal disintegration. The results also indicate that the damage pattern associated with zonal disintegration is influenced by factors including the random distribution type, amplitude of strength fluctuations, the release function, and the release rate due to the excavation. Moreover, the extent of zonal disintegration is governed by in-situ stress conditions and support measures.</p>