## 4 Discussion  
### Covariates Most Relevant to Soil PSF Mapping  
The covariate collection relies on a priori knowledge so that the explicit covariate exhaustive sets are often unlike [21; 26]. Their selection determines the soil mapping performance to a great extent, especially when the number of soil samples is limited, but soil spatial heterogeneity is high [25]. The strategy of this study was to collect environmental covariates as comprehensively as possible, and the selection of key covariates was achieved by covariates reduction techniques.  
AR (all-relevant covariate) sets were used to reveal the underlying process of soil-environment systems of interest [21]. Although the specific covariates retained were differential among AR sets in this study, the Boruta technique retained vegetation and LST (land surface temperature) covariates in both BQ and WQ regions. This indicated that vegetation and LST covariates are the key covariates to identify the spatial variance of soil PSFs of the surface layer. LST is the crucial indicator for the distribution of permafrost [37], and vegetation has an interactive influence on the active layer [60]. Pernafrost degradation on the QTP is in company with the shift in vegetation types and species composition [61; 62], decrease in vegetation biomass, productivity, and species abundance [60; 63; 64]. Under a warming climate, degrading permafrost has profoundly and extensively affected alpine ecology [60]. With the degradation of alpine ecosystem, the soil nutrient declined and surface soil texture became rough [30; 64].  
With different search techniques, minimal-optimal sets were composed of different covariates. It is because the covariates selected by local optimization methods (i.e., GB, GF, and HC) depend on the initial set that the algorithm starts searching with. In general, each data reduction technique has its advantages. For instance, the GF models are the simplest with the least variable redundancy, while GB models retain the most variable information, reflecting more processes [21]. However, they all have their disadvantages as well. The GF models only reflect the major processes and often fell short in model performance. The SA algorithm escapes the local optima by accepting some models that have been degraded due to including or excluding additional sets of variables [52; 65]. In this study, the selected covariates in minimal-optimal sets were different for modeling ALR- and ILR-transformed PSFs. However, they all involved all categories of relevant covariates. For example, in the BQ region, topographic factors were relevant for ALR- and ILR-transformed PSFs. None was retained to predict clay\({}_{\text{irr}}\), while one topographic factor was retained to predict siltl\({}_{\text{irr}}\). Topographic factors were also selected for clay\({}_{\text{irr}}\) prediction while excluded for siltl\({}_{\text{irr}}\) prediction.  
Former studies also emphasized the covariate selection in digital soil mapping, indicating covariate selection can help improve model performance [25; 26]. However, we cannot explicitly assign a universal covariate reduction method for soil PSFs prediction because the predictions based on selected covariates can differ in PSF transformations, modeling techniques, and study areas. It is also the reason why we did not inverse the prediction of transformed PSFs on the basis of the same selection strategy. We prefer to apply and compare multiple covariate selection strategies for each transformed soil fraction based on our results.