The interlayer shearing and crushing fracture zone is one of the key factors on the stability of surrounding rock of underground caverns in interbed layered rock mass. This study focuses on the interlayer shearing and crushing failure of interbed layered rock mass. Interface element was introduced to simulate their contact characteristics. Based on this, the contact model algorithm was improved in this study. In view of the possible contact states for nodes on the boundary of interface, such as bonding contact, slipping contact, opening and embedding contact, the corresponding contact state criterion and modified iterative algorithm were established. It can effectively simulate the nonlinear slip damage between layers of interbed layered rock mass. Then, the proposed model was applied to the stability calculation of an underground cavern of AZAD PATTAN hydropower station in Pakistan. Compared with layered anisotropic constitutive model, the results of the proposed model indicate that the interlayer of interbedded rock mass is prone to mutual dislocation under the action of tectonic stress and excavation load. It is more obvious in the middle and lower parts of the upstream and downstream side walls. The interlayer contact further aggravates the damage range and degree of surrounding rock, and extends to the depth of surrounding rock. The interlayer shear slip failure has an important influence on the stability of high side wall of the cavern. Therefore, the proposed model can reasonably reflect the nonlinear contact characteristics between layers of interbed layered rock mass under tectonic stress. Therefore, it provides a new numerical simulation and analysis method for the stability analysis of underground caverns in interbed layered rock mass.