This article proposes a new methodological approach using a combination of expert knowledge-based (analytic hierarchy process, AHP), bivariate (statistical index, SI) and multivariate (linear discriminant analysis, LDA) models for landslide susceptibility mapping (LSM) in Mazandran Province, Iran. Tolerance and variance inflation factor indicators were used for assessing multi-collinearity among parameters, and three (i.e. profile curvature, soil type and topography wetness index) of 18 factors were eliminated because of multi-collinearity issues. Fifteen geo-environmental conditioning factors including elevation, slope degree, slope aspect, plan curvature, slope length, convergence index, stream power index, distance from river, drainage density, distance from road, distance from fault, lithology, rainfall, land use/landcover and normalized difference vegetation index and 321 landslide locations (testing data set, 70% of total landslides) were used for modeling. The importance of factors showed that distance to road (AHP = 0.201, LDA = 0.301) was the most important factor in landslide occurrence. The validation results using validation data set (138 landslide locations, 30% of total landslides) and area under the receiver operating characteristic curve (AUROC) showed that the ensemble models AHP-LDR (83%), AHP-SI (95%) and SI-LDR (83%) had higher prediction accuracies than the individual AHP (82%), SI (82%) and LDA (79%) models and combination of AHP and SI models along with ALOS-PALSAR remote sensing data and geographic information system (GIS) technique provide powerful tool in LSM in the study area. The results of proposed novel methodological framework can be used by decision-makers and forest engineers for forest management spatially forest roads conservation that have key importance in sustainable development in local and regional scales.