Frequent landslides in the northwestern Himalaya India (NHI) region cause significant loss of life and property, making landslide susceptibility zonation (LSZ) crucial for identifying vulnerable areas. This study aims to develop LSZ maps for the NHI using four bivariate geostatistical models: Frequency Ratio (FR), Weight of Evidence (WoE), Information Value (IV), and Yule's Coefficient (Yc). A total of 38,697 landslides, covering 149.50 km2, were analyzed. The data was split into 70% for training and 30% for testing, ensuring robust model validation. Twelve causative factors were considered, including slope angle, slope aspect, slope curvature, relative relief, terrain roughness index, geomorphon, distance to drainage, land use land cover, lithology, distance to fault/thrust, earthquakes, and rainfall. The models identified high to very high susceptibility zones, covering 28.7%, 32.8%, 48.1%, and 48.2% of the region for the Yc, FR, WoE, and IV models, respectively. ROC analysis revealed that the FR model achieved the highest accuracy, with 0.845 (84.5%) for both validation and prediction. The IV model followed with ROC values of 0.833 (83.3%), while the Yc model performed similarly, with values of 0.831 (83.1%). The WoE model exhibited slightly lower accuracy, with ROC values of 0.830 (83.0%) for validation and 0.831 (83.1%) for prediction. Both the WoE and IV models covered over 98% of landslide areas in high and very high susceptibility zones, indicating a tendency to overestimate highly susceptible areas. The study suggests that the FR and Yc models are particularly effective for LSZ and risk assessment. These results provide valuable insights for hazard management, aiding researchers, planners, and policymakers in selecting appropriate models for LSZ and mitigating landslide risks in other vulnerable regions.