Identifying the spatial variability of soil quality indices is necessary to manage soil resources on a watershed scale. This study aimed to identify suitable indices for soil quality assessment at the watershed scale using various soil characteristics and modeling approaches. Another objective was to map soil quality variability in a representative area in the Qarasu watershed in Kermanshah province, west of Iran. Latin hypercube sampling method using the auxiliary variables used to select 163 sampling points based on land use, soil, and topographical variability in an area of about 57 thousand hectares. In the field operations, soil profiles were described, and samples were taken from different soil profile horizons. Soil properties such as texture, pH, salinity, available water, equivalent calcium carbonate, saturation percentage, soil organic carbon, nitrogen, available phosphorous, potassium, Fe, Zn, Cu and Mn, and soil aggregate stability (mean weight diameter (MWD), geometric mean diametric (GMD), and stable aggregates larger than 0.25 mm (WAS)) measured in the laboratory. Soil quality indices (productivity index (PI), soil quality index (SQI) and reduced dimension soil quality index using principal component analysis (SQI-PCA)) were calculated for each point using the measured soil properties. Soil quality indices were simulated using modeling with the random forest and support vector machine methods and auxiliary variables. Results showed that the range of soil characteristics and integrated soil quality indices was very high across the study area. Soil organic carbon percent varied from about 0.19 to 8.44%. The range of changes in PI in the study area was more than SQI and SQI-PCA indices. Quantities of all soil quality indices were higher in forest and rangeland than in agricultural lands. The spatial estimation accuracy of the SQI-PCA was higher than other soil quality indices and converged well with land use changes compared to PI and SQI.