Abstract: The risk of soil erosion is ever increasing in the Tibetan Plateau, due to the ecological environment under climate warming and human activities. Among them, gully erosion has been the most severe type of soil erosion. Taking the Lhasa River Basin of China as a sample, this study aims to detect the gully occurrence risk using multivariate nonlinear spatial modeling. 2171 gully points were selected to determine the controlling factors of gully erosion using the field survey and remote sensing interpretation. Gully points were set as the value of 1, whereas, no gully point was the value of 0. The Categorical Regression (CATREG), Geodetector, and their combination were utilized to quantify and classify the importance of 15 factors for the risk map of gully occurrence. The 15 factors were the elevation, slope, aspect, topographic position index, topographic wetness index, topographic relief, Normalized Difference Vegetation Index (NDVI), lithology, soil type, permafrost thermal stability, land use, human footprint, distance to the residential area, distance to the river and mean annual precipitation, all of which passed the collinearity test before modeling. Continuous numerical variables were also converted to the type variables by means of discretization. Eight classes were divided by the Jenks Natural Breaks Classification, except for nine classes of aspect. All key factors were converted to the 12.5 m resolution raster datasets. The factor datasets of all sample points were extracted by multi-value extraction to point tool of ArcGIS 10.2 software for data analysis. The Geodetector was run in the Excel software. The factor and risk detectors were used as two parameters of the binary statistical model. The CATREG was performed on the SPSS 20 software. The fitting parameters were obtained, including the factor coefficient, and each class of factors. The improved model was executed in the Raster Calculator Tool in ArcGIS to obtain the map of gully occurrence risk. Finally, the goodness of fit of each was evaluated by the receiver operating characteristic (ROC) curve. The results showed that: 1) The top three factor coefficients in the CATREG were the elevation (0.442), soil type (0.168), and NDVI (0.156). By contrast, the elevation (0.263), soil type (0.251), and human footprint (0.174) were the top three among the Geodetector. 2) The regression performed the best (AUC=0.899), followed by the two combined methods (AUC=0.866/0.848). But, the Geodetector performed relatively poorly (AUC=0.833). All these methods were suitable for spatial modeling in this case. Only a relatively few influences were found in the classification number of factors on the modeling performance. Correspondingly, the classification number with the highest q value of each factor in the Geodetector can be explored to improve the modeling accuracy in the future. In addition, the quantitative score or frequency of adjacent grade can be considered for the grade of sample deficiency. 3) About 9.52% to 13.97% of the Lhasa River Basin was at a very high risk of gully occurrence, mainly distributed in the lower Lhasa River valley and Damxung Basin. 91% of the very high risk areas were clustered in the elevation classes 1-3, whereas, 79% of the high risk areas were clustered in the elevation classes 2-5. Therefore, the hydrological connectivity can be altered to restore the vegetation for less topsoil disturbance in the low-elevation area around large towns. The high-elevation area with a low risk can be used to prevent overgrazing and vegetation destruction because of the high precipitation in the northeast of the basin. The finding can provide a strong reference to construct the ecological security barrier, as well as the soil and water conservation on the Tibetan Plateau.