<p>Loess-mudstone interface landslides occur frequently on the Chinese Loess Plateau, where rainfall serves as the primary natural trigger of such failures. This study aims to investigate the seepage and stability characteristics of this type of slope under rainfall conditions. Complex preferential seepage channels commonly develop within such slopes and can be classified into sequential, non-sequential, and interface types. To this end, a bidirectionally coupled modeling framework integrating surface runoff and subsurface seepage is proposed. This framework enables a coordinated representation of the hydraulic behaviors of various preferential seepage channels, accurately capturing the seepage characteristics of slopes. The hydrological module is sequentially coupled with a soil mechanics model to analyze slope stability. The model is numerically implemented using the Partial Differential Equation (PDE) module of COMSOL Multiphysics. The model is validated against a laboratory experiment and subsequently applied to a representative slope in Gansu Province, China. The results show that under short-duration, high-intensity rainfall, non-sequential preferential flow develops significantly, forming a high perched water level and correspondingly high positive pore water pressure above the low-permeability interface. Combined with stress concentration at the crack tip, the rear-edge tensile crack tends to propagate and coalesce with the loess-mudstone interface. By contrast, under long-duration, low-intensity rainfall, interface preferential flow dominates, making the slope more susceptible to sliding failure along the interface. The proposed framework provides an effective approach for analyzing the seepage characteristics and stability evolution of most types of natural slopes under rainfall conditions.</p>