<p>Slopes are considered porous media in rainfall infiltration analysis based on two-phase flow model, while the impact of preferential flow through internal fractures is typically overlooked. This study aims to investigate water and air transport patterns within fractures during rainfall and their impact on the stability evolution. In this paper, a novel framework is developed to couple two-phase flow with the stability analysis of fractured rock slopes. Firstly, a hydromechanical coupled model is proposed and verified to simulate water‒air interactions. Secondly, two original FISH functions are proposed to implement a more physically representative relative permeability model and a more realistic boundary condition for simulating rainfall infiltration under two-phase flow model. Further, the failure modes with different fracture distribution characteristics are revealed using the Universal Distinct Element Code. Finally, the applicability of the proposed model is demonstrated by a practical case study at the GS Hydropower Station. The results indicate that air entrapment within fractures significantly delays rainwater infiltration, accelerating the transition from the flow boundary to the pressure boundary at the infiltration end of the fractures. This ultimately mitigates the progression of slope instability. Additionally, the choice of relative permeability models is found to critically influence both rainfall infiltration behaviour and the subsequent slope stability evolution. The results may provide a helpful reference for hazard assessment and control of rainfall-induced landslides in fractured rock slopes.</p>