Assessing Variation of Soil Quality in Agroecosystem in an Arid Environment Using Digital Soil Mapping

Abstract

Monitoring the soil quality (SQ) in agricultural ecosystems is necessary for using sustainable soil and land resources. Therefore, to evaluate the SQ variation in an arid environment in the Bajestan region, northeastern Iran, two soil quality indices (weighted additive soil quality index-SQI_w and nemoro soil quality index-SQI_n) were applied. SQIs were assessed in two datasets (total data set-TDS and minimum data set-MDS) by linear (L) and nonlinear (NL) scoring methods. Physicochemical properties of 223 surface soil samples (0–30 cm depth) were determined. The random forest (RF) model was used to predict the spatial variation of SQIs. The results showed the maximum values of the SQIs in areas with saffron land covers, while the minimum values were acquired in the north of the study area where pistachio orchards are located due to higher EC and SAR. The environmental variables such as topographic attributes and groundwater quality parameters were the main driving factors that control SQIs distribution. These findings are beneficial for identifying suitable locations sites to plan agricultural management and sustainable usage of groundwater resources strategy to avoid further increase of soil salinity.

1. Introduction

Due to the limitation of arable land in the arid and semiarid area and the lack of water, increasing production by adding the area under cultivation is impossible [1]. Therefore, meeting the nutritional needs of this growing population requires optimal use of agricultural land and water, where inappropriate land use and soil management cause land degradation [2]. Soil is recognized as a critical component of ecosystem sustainability, which is required for long-term development and the most efficient use of natural resources [3]. On the other hand, proper identification of different soil physicochemical properties has a significant role in determining the erodibility, soil degradation, and management of agricultural lands and protection of soils, especially in arid and semiarid regions [4,5,6,7]. Soil quality (SQ) assessment is a practical approach to identifying the primary effects of management practices [8]. Understanding SQ is important for identifying problem areas and evaluating sustainable agricultural management [1,4].

The SQ of the agricultural lands is affected by different environmental and management factors, including landform and topographic situation [2,6], land use/cover [5], irrigation water quality [2,6,9,10], the amounts of inputs and outputs of soil organic carbon (SOC) and chemical fertilizers [11,12,13]. Several scholars [6,7] stated that intensive tillage operations with the removal of crop residues, excessive land use, and failure to meet the need for soil fertilizers, the use of poor quality irrigation water could lead to reduced production and SQ in arid and semiarid regions. Some soil properties (e.g., SOC, available phosphorous and potassium, soil salinity, and alkalinity) were commonly applied as an appropriate indicator for identifying SQ in Iran [3,8,11,14,15,16].

So far, various methods have been used to collect data, and measure and evaluate SQ from mostly qualitative to quite quantitative methods, including scorecards, visual soil assessments, field kits, geostatistical methods and laboratory analysis [17]. Parameters affecting soil quality indicators (SQI) in the form of soil processes and properties are defined as being sensitive to changes in soil management. These properties can be a set of physical, chemical, biological or a combination of these properties. Different sets of characteristics were proposed to affect the SQ and determine the SQI [17,18]. Scholars have determined the SQI based on the total data set (TDS) and minimum data set (MDS) of properties affecting SQ [3,7,8,11,13,14,15,16].

Because it is complicated to interpret many variables and draw conclusions from them, it is recommended that the set of variables be combined into one index. This is done by combining the data and applying the appropriate weight to each variable. Currently, many quantitative indicators such as weighted additive SQI (SQI_w) and Nemro SQI (SQI_n) have been developed to calculate SQI [8,11,16].

Many studies of SQI assessment have been done in different parts of the world; however, a few studies have been done to study the effect of different environmental covariates (i.e., topographic attributes, groundwater quality parameters, remote sensing indices) through agricultural lands on SQI using digital soil mapping (DSM, [2,11,16,19]. The “scorpan” model [20] was applied to predict the spatial variability of soil properties based on the relationship between environmental covariates and soil properties for DSM studies [21,22,23]; then, a significant number of linear and nonlinear techniques were applied to predict object property.

Agricultural lands in eastern Iran on the playas margin suffer from soil salinity and alkalinity, but, by improvement practices, the area under cultivation here has increased. In recent years, the cultivation land of pistachio has increased in the Bajestan region, especially in the playa margin, in which some parts of this area showed a reduction in fertile lands and irrigation water quality [6]. Determining the quality of these lands to continue their use and determining SQI for areas rehabilitated by cultivation can be important and influential for managing long-term agricultural use of lands within playa margins without reduction SQ. The objectives of the present study were to (i) explore the SQ variation through pistachio, pomegranate, saffron, and barley in the cultivation lands; (ii) provide a framework to assess SQ using different SQI scoring functions (i.e., linear and nonlinear) and DSM technique for future land use planning; and (iii) spatially predict SQIs maps and determine the importance environmental covariates affecting SQ for a better understanding of the land potential and suitability.

2. Materials and Methods

2.1. Study Area

The study was conducted in the Bajestan region, southwest of Khorasan Razavi (57°57′56″ to 58°00′40″ E, and 34°17′91″ to 34°33′79″ N, Figure 1) with an area of about 172,419 ha. It is characterized by an arid climate with mild winters and dry, hot summers. The mean annual temperature and precipitation are 17.3 °C and 193 mm, respectively. The elevation ranges from 786 to 2283 m above sea level (a.s.l.) in the studied area. This region has Aridic moisture and Thermic temperature regimes [24]. Rainfed agriculture is predominant in this region with major crops including pomegranate (Punica granatum L.), pistachio (Pistacia vera L.), barley (Hordeum vulgare L.) and saffron (Crocus sativus L.) which are located on different landforms, including pediments, alluvial fans, playa clay flats and dune fields (Figure 1b). In general, pistachio orchards and barley cultivated lands are found in the north part of the study areas located in playa margins, while pomegranate and saffron are located in the center and south of the area (Figure 1b). Overall, the area’s geologic material through the mountain to flatlands consists of lower Cretaceous undifferentiated rocks, partly massive and bedded limestone, and recent and old alluvial deposits [6]. According to previous studies in the study area, most soil orders in the study area are Typic Haplocalcids (WRB: Haplic Calcisols (Loamic, Raptic) over Luvic Skeletic Calcisols (Arenic, Raptic)), Typic Torrifluvents (WRB: Eutric Fluvisols (Loamic) over Luvic Skeletic Calcisols (Arenic, Raptic)) and Typic Torripsamments (WRB: Calcaric, Eutric Arenosols (Aeolic)) [25].

2.2. Field Sampling and Soil Properties Analysis

A total of 223 soil samples were collected from the 0–30 cm depth. It should be noted that each soil sample is a combination of 5 soil subsamples around the main locations. All soil samples were passed through a 2 mm sieve after air-drying. The titration method was used to determine calcium carbonate equivalent (CCE) [26]. SOC was measured by Walkley–Black method [27]. The soil pH and electrical conductivity (EC) were determined in saturated soil paste and extracted of saturated soil paste using pH meter [28] and EC meter [29], respectively. The saturation percentage (SP) was measured in soil paste [29]. The soil texture was determined based on the hydrometer method [30]. Total nitrogen (TN) was measured by the Kjeldahl approach [31], and the available phosphorus (P_av) was determined by the colorimetric method [32]. Available potassium (K_av) was extracted with 1 M NH₄OAC (pH 7) and measured via flame photometer [33]. Soluble calcium (Ca_aq) and magnesium (Mg_aq) were measured by titration method with EDTA and soluble sodium (Na_aq) by flame photometry [34]. Sodium adsorption ratio (SAR) was calculated by the following equation [35]:

$$SAR=\frac{N{a}_{aq}}{\sqrt{\frac{C{a}_{aq}+M{g}_{aq}}{2}}}$$

(1)

2.3. Soil Quality Index (SQI) Assessment

2.3.1. Total Data Set (TDS) and Minimum Data Set (MDS)

The TDS and MDS were used to calculate the SQI. In TDS, all available soil properties were generally used to provide a comprehensive overview of SQ in the study area. However, the MDS was used to reduce the volume of soil property data, and the key indicators were selected [11,16,17]. In this study, all accessible soil properties (i.e., EC, pH, SP, CCE, SOC, TN, P_av, K_av, Ca_aq, Mg_aq, Na_aq, and SAR) were included in the TDS method. The selected soil properties define soil health, productivity, fertility, soil degradation and soil and water interaction. Principal component analysis (PCA) was applied to reduce dimensionality in the data set and determine the most important properties to include in the MDS [18]. The factors with eigenvalues ≥1 and soil properties with the highest loadings in each PC were presumed as the best indicator of SQ. However, when more than one soil variable with the highest loadings is retained within a PC, Pearson’s correlation coefficients among the soil properties were used to detect any redundant variable. In the MDS, the feature with the highest value is selected among the features that significantly correlate with each other. For attributes selection in the MDS, between highly-correlated attributes (r > 0.7), the variables with the highest factor loading were considered [18].

2.3.2. Indicator Scoring

Linear (L) and nonlinear (NL) methods were used to transform soil properties into a dimensionless score (between 0.1 and 1) using the following functions: ‘more is better’ for soil properties such as SP, Ca_aq, Mg_aq, SOC, TN, P_av, K_av; ‘less is better’ for SAR, Na_aq, CCE; and optimal range for pH and EC to the TDS and MDS [5,16] (Table S1). The optimal values 0.2–2 dS m⁻¹ and 7 were used for EC and pH, respectively [36,37].

2.3.3. Weight Assignment and SQIs

We used two SQI equations, including nemoro SQI (SQI_n) and weighted additive SQI (SQI_w) [5,14,15] as following equations:

$$SQ{I}_n=\sqrt{\frac{P_{\mathit{ave}}^{}{}^2+P_{\min }^{}{}^2}{2}}\times \frac{n-1}{n}$$

(2)

$$SQ{I}_w={\displaystyle \sum _{i=1}^n{W}_i}{N}_i$$

(3)

where N_i, n, P_ave, P_min, and W_i are the indicator scores, the number of indicators, the average value, the minimum value for the scores attained for each sampling point, and the weight assigned to each indicator, respectively. The weights in the SQI_w approach were obtained according to the communality of PCA in both TDS and MDS. The weights of every indicator were determined based on the ratio of the communality [38].

2.3.4. Soil Quality (SQ) Grades and Comparison of Indices

In this study, eight different SQIs were calculated based on different data sets (TDS and MDS), SQI approaches (SQI_w and SQI_n), and two scoring systems (L and NL). Based on Johnson and Wichern’s [39] definition, each SQI was classified into five classes, including very high (I), high (II), moderate (III), low (IV) and very low (V). The performance of different SQIs was evaluated by match/mismatch analysis. Thus, for assessment of the agreement level between the soil grades (e.g., very high, high, moderate, low and very low), the Kappa statistic was applied [8,16]. Therefore, based on Nabiollahi et al. [8] and Zeraatpisheh et al. [16], the Kappa statistic agreement levels were classified into six levels, including none (<0), poor (0–0.19), weak (0.20–0.39), moderate (0.40–0.59), strong (0.60–0.79) and excellent (0.80–1.00). Additionally, the correlation analysis was performed between the indices to better understand the relationships among the SQIs. The relationship for indicator methods was examined by regression analysis using the SAS software version 9.4 [40].

2.4. Spatial Prediction of SQIs

In order to predict and map the spatial distribution of SQIs in the study area, the random forest (RF) model through DSM was used to assess the relationship between SQIs and environmental covariates [8,16,20]. RF was considered one of the most accurate, efficient, and popular machine learning algorithms [19,22,41]. A set of environmental covariates explaining the “scorpan” factors [20] were used as predictors in the study area to predict SQI distributions (Table S2). A total of 21 terrain attributes, including the first and second digital elevation model (DEM) derivatives (28 m × 28 m), were extracted from an SRTM 1 arc-second DEM (http://earthexplorer.usgs.gov, accessed on 15 July 2020). The SAGA GIS software (System for Automated Geoscientific Analysis) version 2.2 was used to obtain terrain covariates [42]. Moreover, the remotely sensed data and 52 RS indices were obtained from Sentinel-2 (https://sentinel.esa.int/web/sentinel/missions, accessed on 15 July 2020), including high-resolution spatial and temporal information with 13 bands with spatial resolutions of 10 m, 20 m and 60 m [43]. Additionally, the groundwater quality parameters (e.g., HCO₃⁻, Cl⁻, SO₄²⁻, Na⁺, Ca²⁺, Mg²⁺, EC, pH, SAR and TDS), geology, geomorphology and land use maps were used as predictive variables (Table S2). A total of 190 groundwater resources, including Qanat and wells, were selected and sampled from those water for analysis of groundwater quality parameters according to standard procedures [44].

Regarding many environmental covariates, the RF model provides variable importance ranging from 0 to 100%. When the relative importance of the environmental covariates was below 15%, they were considered unimportant covariates and were removed from the models [45]. Finally, the models were trained by the rest of the environmental covariates. Therefore, for modeling different SQIs, 25–36 environmental covariates were used in the random forest model.

The spatial modeling of the SQI was done using the RF in the “caret” package in R3.3.1 [46]. The accuracy of the SQ value maps was validated using three validation criteria: the root mean square error (RMSE); the mean absolute error (MAE); and the coefficient of determination (R²) using 10-times leave-one-out cross-validation [47]. The Kappa statistic was used to examine the accuracy of the soil grade classification by the RF model compared to the observed SQ grades [8,16].

2.5. Variable Importance for Soil Quality Indicators (SQIs) Maps

Variable importance helps select important factors in predicting and modeling soil properties [48]. The “Out-of-Bag (OOB)” error (Equation (4)) is estimated randomly by changing the values of the covariates, which is considered as variable importance. The variable importance varied between 0–100%.

$$\mathrm{OOB} \mathrm{error}=\frac{1}{N}{\displaystyle \sum _{i=1}^{N}I\left[Y_{OOB}\left(X_i\right) \ne Y_i\right]}$$

(4)

where I [Y_OOB(X_i) ≠ Y_i] is recognized as an indicator function equal to 0 when the predicted and actual classes are the same and equal to 1 otherwise.

According to the previous studies, the participation rate of covariates and their importance were determined using the RF model in the “caret” package in R3.2.5 [46], which produces a reliable spatial prediction of soil properties [21,47]. However, the importance of a variable depends on the applied method. All input environmental factors are introduced as model random variables, representing their uncertainty. This test can be used to access contributions by the factors (i.e., relative importance) to the model.

2.6. Statistical Analysis

The mean values of SQIs and individual soil properties were used to define the significant differences (p < 0.05) among the land covers using the analysis of variance (ANOVA) approach in SAS software [40].

3. Results and Discussion

3.1. Descriptive Statistics of Soil Properties

The summary descriptive statistics of soil properties used to assess SQIs in this study are presented in

Table 1. Descriptive statistics of soil properties in the study area (n = 223).

Variable	Unit	Mean	Minimum	Median	Maximum	StDev	CV%	Skewness	Kurtosis
EC	(dS m−1)	7.19	0.14	5.58	46.60	6.11	84.94	2.31	9.39
pH	-	7.85	6.10	7.85	9.62	0.51	6.44	−0.13	1.99
SP	(%)	30.35	18.00	29.10	54.23	5.83	19.21	0.66	0.94
CCE	(%)	17.55	2.00	17.25	45.00	7.59	43.26	1.10	3.06
SOC	(%)	0.59	0.02	0.35	3.06	0.58	99.24	1.54	2.04
TN	(%)	0.06	0.00	0.03	0.91	0.08	135.39	6.62	68.86
Pav	(mg kg−1)	11.64	0.20	8.00	80.00	12.06	103.66	2.43	7.92
Kav	(mg kg−1)	215.40	17.00	199.00	695.00	107.68	49.99	1.16	1.96
Caaq	(meq L−1)	21.91	1.50	21.40	73.60	13.54	61.80	0.83	1.05
Mgaq	(meq L−1)	10.54	0.60	8.00	56.00	7.94	75.31	2.50	9.16
Naaq	(meq L−1)	41.65	0.46	26.80	272.00	43.03	103.33	1.68	4.00
SAR	-	11.06	0.46	9.00	41.00	9.44	85.32	0.86	0.04

EC—electrical conductivity; pH—soil reaction; SP—saturation percentage; CCE—calcium carbonate equivalent; SOC—soil organic carbon; TN—total nitrogen; P_av—available phosphorous; K_av—available potassium; Ca_aq—soluble calcium; Mg_aq—soluble manganese; Na_aq—soluble sodium; SAR—sodium absorption ratio.

. The results demonstrated that the range of soil properties were: 0.14 to 46.60 (dS m⁻¹) for EC; 6.10 to 9.62 for pH; 18.0 to 54.23 (%) for SP; 2.0 to 45.0 (%) for CCE; 0.02 to 3.06 (%) for SOC; 0 to 0.91 (%) for TN; 0.20 to 80.0 (mg kg−1) for P_av; 17.0 to 695.0 (mg kg⁻¹) for K_av; 1.50 to 73.60 (meq L⁻¹) for Ca_aq; 0.60 to 56.0 (meq L⁻¹) for Mg_aq; 0.46 to 272.0 (meq L⁻¹) for Na_aq; and 0.46 to 41.0 for SAR (Table 1). The lowest and highest coefficient of variation (CV) was also observed in pH and TN, respectively (Table 1). According to the CV classification of Pahlavan-Rad and Akbarimoghaddam [48], pH and SP indicated low and medium CV classes, respectively, while other soil properties had high and very high variability (Table 1).

As shown in Figure 2, the highest positive correlations were observed between Na_aq and EC (r = 0.78), SAR and EC (r = 0.67), SAR and Na_aq (r = 0.85), Mg_aq and Ca_aq (r = 0.58), Mg_aq and Na_aq (r = 0.47), TN and SOC (r = 0.63), P_av and SOC (r = 0.54) while the highest significant negative correlations were found between pH and CCE (r = −0.28), SOC and EC (r = −0.27), SOC and Na_aq (r = −0.25), SOC and SAR (r = −0.26).

3.2. Variation Changes of SQ through Different Cultivated Lands

3.2.1. TDS Indicator Method

Table 2. Communality and weight of SQIs in the TDS and MDS methods.

Indicator	TDS		MDS
	COM a	Weight	COM	Weight
EC (dS m−1)	0.785	0.079	-	-
pH	0.893	0.090	0.676	0.128
SP (%)	0.773	0.078	-	-
CCE (%)	0.857	0.086	0.875	0.166
SOC (%)	0.828	0.084	0.821	0.156
TN (%)	0.799	0.081	-	-
Pav (mg kg−1)	0.691	0.070	-	-
Kav (mg kg−1)	0.774	0.078	-	-
Caaq (meq L−1)	0.859	0.087	0.990	0.188
Mgaq (meq L−1)	0.797	0.080	-	-
Naaq (meq L−1)	0.929	0.094	0.940	0.179
SAR	0.921	0.093	0.959	0.182

^a Communality. EC—electrical conductivity; pH—soil reaction; SP—saturation percentage; CCE—calcium carbonate equivalent; SOC—soil organic carbon; TN—total nitrogen; P_av—available phosphorous; K_av—available potassium; Ca_aq—soluble calcium; Mg_aq—soluble manganese; Na_aq—soluble sodium; SAR—sodium absorption ratio; TDS—total data set; and MDS—minimum data set.

is presented the PCA results. The results of PCA illustrated which Na_aq (weight = 0.094) and P_av (weight = 0.070) had the highest and lowest weight among all studied soil properties, respectively when the PCA was applied based on TDS (Table 2).

Table 3. Comparison of the mean values of soil properties through different cultivated lands.

Land Cover	Pistachio	Barley	Pomegranate	Saffron	Pr > F
N	99	29	50	45	-
EC (dS m−1)	9.58 a	9.44 a	4.13 b	3.88 b	0.0001 **
pH	7.74 b	7.87 ab	7.97 a	7.95 ab	0.0283 *
SP (%)	29.33 b	27.35 b	32.60 a	32.02 a	0.0001 **
CCE (%)	17.69 ab	20.04 a	15.82 b	17.54 ab	0.0237 *
SOC (%)	0.29 c	0.45 c	1.16 a	0.70 b	0.0001 **
TN (%)	0.03 c	0.04 bc	0.11 a	0.06 b	0.0001 **
Pav (mg kg−1)	7.16 c	5.97 c	21.65 a	14.00 b	0.0001 **
Kav (mg kg−1)	197.10 cb	165.00 c	260.80 a	237.70 ab	0.0001 **
Caaq (meq L−1)	23.63 a	22.59 a	16.14 b	24.10 a	0.007 **
Mgaq (meq L−1)	11.21 a	10.27 a	10.80 a	8.94 a	0.4577 ns
Naaq (meq L−1)	59.52 a	51.66 a	22.91 b	16.69 b	0.0001 **
SAR	15.19 a	13.42 a	6.72 b	5.27 b	0.0001 **

** and * denote significance at the 0.01 and 0.05 probability levels, respectively, while ns is insignificant. Means with the same letter are not significantly different based on Duncan’s multiple range test. EC—electrical conductivity; pH—soil reaction; SP—saturation percentage; CCE—calcium carbonate equivalent; SOC—soil organic carbon; TN—total nitrogen; P_av—available phosphorous; K_av—available potassium; Ca_aq—soluble calcium; Mg_aq—soluble manganese; Na_aq—soluble sodium; SAR—sodium absorption ratio.

and

Table 4. Comparison of the mean values of soil quality indices through different cultivated lands.

Land Cover	Pistachio	Barley	Pomegranate	Saffron
N	99	29	50	45
SQI-w-L-TDS	0.456 b	0.455 b	0.539 a	0.523 a
SQI-w-NL-TDS	0.369 b	0.370 b	0.562 a	0.539 a
SQI-n-L-TDS	0.300 c	0.300 c	0.361 a	0.345 b
SQI-n-NL-TDS	0.238 b	0.240 b	0.370 a	0.355 a
SQI-w-L-MDS	0.544 b	0.555 b	0.626 a	0.628 a
SQI-w-NL-MDS	0.374 b	0.386 b	0.529 a	0.564 a
SQI-n-L-MDS	0.333 b	0.374 b	0.391 a	0.392 a
SQI-n-NL-MDS	0.223 b	0.238 b	0.328 a	0.339 a

Means with the same letter are not significantly different based on Duncan’s multiple range test (significant at the 0.01 probability levels). SQI_w—weighted additive soil quality index; SQI_n—nemoro soil quality index; TDS—total data set; MDS—minimum data set; L—linear; and NL—non-linear.

show the mean comparison of the soil properties and SQIs in different land covers in the studied area. Among all soil properties studied, Mg_aq did not show significant differences between different land covers (Table 3). The highest amount of EC (9.58 dS m⁻¹), Ca_aq (23.63 meq L⁻¹), Mg_aq (23.63 meq L⁻¹), Na_aq (23.63 meq L⁻¹), and SAR (15.19) were found under cultivated lands with pistachio orchards while soils under pomegranate orchards includes of the highest values of pH (7.97), SP (32.60%), SOC (1.16%), TN (0.11%), P_av (21.65 mg kg⁻¹), and K_av (260.80 mg kg⁻¹). These results are in line with Maleki et al. [6], who reported the increasing values of base cations in soils under pistachio orchards through 14 years in Bajestan cultivated lands. The maximum value of CCE was found in soils under barley cultivated lands (Table 3). However, there is no significant difference in importance between different land covers. According to an earlier study in this area by Maleki et al. [6], most cultivated lands are located on the low calcareous parent material.

It should be noted that the mean value of all SQIs calculated based on L scoring had better results than NL scoring (Table 4). The results show that pomegranate and saffron cultivated lands have the highest amount of SQI_w and SQI_n in the TDS data set, which directly results in higher values in P_av, K_av, SOC, TN, SP and lower values in EC, SAR, Na_aq in the mentioned land covers (Table 3). The pistachio and barley cultivated lands have the lowest SQ in the TDS data set and significantly differ from the two land covers, including pomegranate and saffron. The higher values of EC, SAR, and Na_aq are indicated under pistachio and barley cultivated lands that decrease SQ [3,17]. Nabiollahi et al. [8] and Castro et al. [12] also reported the importance of soluble salts concentration and soil salinity on decreasing SQ in arid and semiarid regions.

One of the main soil properties that have a significant effect on increasing SQ is SOC content [12,49], which affects most soil physicochemical (i.e., TN, soil stability, structure, infiltration, and water retention) and biological properties [13,16,50]. Therefore, by taking management measures to increase SOC, other properties can be optimized, increasing the SQ in pistachio and barley cultivated lands. In pistachio orchards and barley lands, low litter amounts and nonreturn of plant residuals to the soil, respectively, reduce the amount of surface cover and the quality and quantity SOC, and consequently, SQ has decreased [10,11,13,51]. In this regard, Nie et al. [51] showed that without adding organic matter and in the control treatment, SQ was at the lowest level and was classified in class 4, but with the addition of organic matter, SQ class was upgraded by one to two degrees.

In the comparison between the two methods, SQI_w and SQI_n, it can be seen that the numerical value of SQI_n in both TDS and MDS is less than the SQI_w (Table 4). Because SQI_w for the studied soil properties and scoring also considers weight, while in the SQI_n, the score is calculated only based on the average values and the minimum score of the properties. Many other researchers have concluded that the value of SQI_n is lower than the SQI_w in all conditions [14,16].

3.2.2. MDS Indicator Method

Table 2 and

Table 5. Results of principal component analysis of eight soil properties in the studied area.

PCs a	PC1	PC2	PC3	PC4	PC5
Eigenvalue	3.46	2.08	1.38	1.23	0.99
Percent	28.90	17.30	11.50	10.20	8.10
Cumulative percent	28.90	46.20	57.70	68.00	76.10
Eigenvectors
EC (dS m−1)	0.43	0.23	0.16	0.03	0.06
pH	−0.14	0.03	0.43	−0.48	−0.06
SP (%)	−0.15	0.31	0.19	0.40	−0.40
CCE (%)	0.15	0.07	−0.40	0.35	0.58
SOC (%)	−0.31	0.40	−0.01	−0.28	0.21
TN (%)	−0.27	0.35	−0.03	−0.27	0.31
Pav (mg kg−1)	−0.28	0.38	−0.04	0.18	0.17
Kav (mg kg−1)	−0.12	0.38	0.18	0.43	−0.21
Caaq (meq L−1)	0.25	0.24	−0.44	−0.18	−0.44
Mgaq (meq L−1)	0.26	0.37	−0.35	−0.28	−0.15
Naaq (meq L−1)	0.45 bc	0.25	0.22	−0.05	0.07
SAR	0.40	0.13	0.44	0.03	0.26

^a Principal Component. ^b Bold factor loadings selected as MDS. ^c Underlined factor loadings are considered highly weighted. EC—electrical conductivity; pH—soil reaction; SP—saturation percentage; CCE—calcium carbonate equivalent; SOC—soil organic carbon; TN—total nitrogen; P_av—available phosphorous; K_av—available potassium; Ca_aq—soluble calcium; Mg_aq—soluble manganese; Na_aq—soluble sodium; SAR—sodium absorption ratio.

showed that the soil properties with the highest factor loading in each PC were selected for MDS. EC, SP, TN, P_av, K_av, and Mg_aq were not used in MDS because they had lower factor loading in each PC than other soil properties (Table 5). The remaining soil properties are applied in the MDS (Table 2). According to Table 5, the first four PCs demonstrated eigenvalues > 1, explaining 68.0% of the total variance, and the fifth PCs around 76.10%. The CCE, SOC, Ca_aq, Na_aq, and SAR showed communalities > 82% of MDS variance (Table 2). Additionally, the communalities for pH in MDS described 67% of the variance through five components. Overall, using the MDS approach could be a reliable method to reduce the number of soil properties, consequently resulting in reduced time consumption, workload, and related costs [8,11,16]. Comparing NL and L scoring system results for MDS, the results showed higher values for the L scoring system (Table 4).

3.3. Assessment of SQ Grades through Different Cultivated Lands

3.3.1. TDS Method

Both SQI_w and SQI_n were classified into five categories using both scoring systems (

Table 6. Soil quality grades classification for indices and indicator methods.

Index	Indicator Method	SSF	Soil Quality Grades
			I (Very High)	II (High)	III (Moderate)	IV (Low)	V (Very Low)
SQIw	TDS	Linear	>0.568	0.506–0.568	0.444–0.506	0.382–0.444	<0.382
	MDS	Linear	>0.691	0.608–0.691	0.525–0.608	0.442–0.525	<0.442
	TDS	Non-linear	>0.639	0.529–0.639	0.418–0.529	0.307–0.418	<0.307
	MDS	Non-linear	>0.718	0.575–0.718	0.432–0.575	0.289–0.432	<0.289
SQIn	TDS	Linear	>0.393	0.347–0.393	0.303–0.347	0.258–0.303	<0.258
	MDS	Linear	>0.468	0.407–0.468	0.347–0.407	0.286–0.347	<0.286
	TDS	Non-linear	>0.447	0.367–0.447	0.288–0.367	0.208–0.288	<0.208
	MDS	Non-linear	>0.478	0.382–0.478	0.287–0.382	0.191–0.287	<0.191

SQI_w—weighted additive soil quality index; SQI_n—nemoro soil quality index; TDS—total data set; MDS—minimum data set; SSF—standard scoring functions.

). The SQI_w and SQI_n maps were prepared based on the TDS and L and NL scoring systems (Figure 3A,B), which indicates that the dominant parts of the region are in high (II) and moderate-quality (III) classes. Soils under pomegranate and saffron cultivated lands showed the highest grades in high (II) classes for SQI_w in both L and NL scoring systems (Figure 3A and

Table 7. Percentage of soil quality grades for different indices, data sets (TDS and MDS), and scoring methods in the studied area.

Index	Data Set	Land Cover	Very High (I) *		High (II)		Moderate (III)		Low (IV)		Very Low (V)
			Linear **	Non-Linear	Linear	Non-Linear	Linear	Non-Linear	Linear	Non-Linear	Linear	Non-Linear
SQIw	TDS	Total area	9.42	10.76	27.80	21.52	41.70	22.42	16.59	22.87	4.48	22.42
		Pomegranate	4.93	5.38	11.66	10.76	5.83	4.48	0.00	1.79	0.00	0.00
		Pistachio	0.00	0.45	6.28	2.69	23.32	10.76	11.21	14.35	3.59	16.14
		Saffron	4.48	4.48	8.07	7.17	6.73	4.93	0.90	2.69	0.00	0.90
		Barley	0.00	0.45	1.79	0.90	5.83	2.24	4.48	4.04	0.90	5.38
	MDS	Total area	6.28	4.04	31.39	21.08	39.01	29.15	18.39	25.11	4.93	20.63
		Pomegranate	1.35	0.45	13.00	7.17	8.07	10.76	0.00	4.04	0.00	0.00
		Pistachio	0.90	0.90	7.17	2.24	18.83	12.11	13.45	13.45	4.04	15.70
		Saffron	3.14	2.24	9.87	9.87	5.38	4.04	1.79	2.24	0.00	1.79
		Barley	0.90	0.45	1.35	1.79	6.73	2.24	3.14	5.38	0.90	3.14
SQIn	TDS	Total area	5.38	4.93	23.32	21.08	37.67	20.18	27.35	29.60	6.28	24.22
		Pomegranate	4.04	1.35	11.66	13.00	5.83	4.48	0.90	3.59	0.00	0.00
		Pistachio	0.00	0.45	1.79	0.45	19.73	8.52	17.94	17.94	4.93	17.04
		Saffron	1.35	3.14	8.52	6.28	8.52	5.38	1.79	4.48	0.00	0.90
		Barley	0.00	0.00	1.35	1.35	3.59	1.79	6.73	3.59	1.35	6.28
	MDS	Total area	3.59	4.04	16.14	8.97	35.87	30.04	35.87	28.70	8.52	28.25
		Pomegranate	2.24	0.45	4.04	3.59	12.11	10.31	4.04	8.07	0.00	0.00
		Pistachio	0.00	0.90	3.59	0.90	13.45	9.42	20.18	12.56	7.17	20.63
		Saffron	0.90	2.24	6.73	3.59	7.62	8.52	4.93	4.04	0.00	1.79
		Barley	0.45	0.45	1.79	0.90	2.69	1.79	6.73	4.04	1.35	5.83

* Grade; ** Scoring method. SQI_w—weighted additive soil quality index; SQI_n—nemoro soil quality index; TDS—total data set; MDS—minimum data set.

). In contrast, barley and pistachio orchards are in the moderate (III) to low (IV) grades (Table 7). Figure 3A,B also shows that most areas with very low (V) SQ are located in the northeastern and northwestern regions, where the largest area of pistachio and barley cultivation areas is recognized (Figure 1). The reason for the low SQ in these areas can be attributed to the higher EC and SAR in areas under pistachio and barley cultivation. High soil salinity and alkalinity values were attributed to primary salinity factors such as proximity to the playa and secondary factors such as low irrigation water quality (for more detailed information, see Maleki et al. [6].

3.3.2. MDS Method

The SQI_w and SQI_n in both L and NL scoring systems showed lower SQ results for MDS than TDS (Table 6). The results of Figure 3A,B, and Table 7 indicated that most of the study areas using the SQI_w (i.e., both L and NL scoring systems) had grades II and III. However, the SQI_n maps (i.e., both of L and NL scoring system) indicate that a large area of the region’s north is in classes IV and V (Figure 3A,B). As previously discussed in TDS, areas under pistachio and barley land covers had the most dominant soil grades III and IV (Figure 3A,B, and Table 7).

3.4. Indices Comparison and Evaluations

The strong linear relationship of SQI between the two sets of TDS and MDS for all soil samples in the study area is shown in Figure 4. As can be seen, R² between MDS and TDS of SQI_w and SQI_n for 0.79 and 0.75, respectively (Figure 4a), by linear scoring while R² ranged from 0.83 to 0.80 for SQI_w and SQI_n in NL scoring (Figure 4b). On the other hand, the results of R² between SQIs using L and NL ranged from 0.77 to 0.84 (Figure 5a,b). These findings illustrated and emphasized its reality and accuracy for SQ that it is possible to confidently use MDS instead of the TDS in different regions as same as studied area. Several studies have reported an excellent explanatory of R² between these SQIs (e.g., [8,13,16]. Therefore, the various studies showed an established relationship between TDS and MDS. While the TDS model is more accurate due to using soil properties, the MDS approach is more economical as it used fewer soil properties.

The results of Pearson’s correlation coefficients among all of the different SQIs exhibited a strong correlation (Figure 6) as ranging from 0.77 to 0.99. The lowest correlation was identified between SQI_w-NL-MDS vs. SQI_n-L-TDS (r = 0.77) and SQI_n-NL-MDS vs. SQI_n-L-TDS (r = 0.77).

Based on the Kappa statistic agreement reported by Emami et al. [52], the Kappa statistical analysis in TDS and MDS showed a weak to moderate satisfactory level of agreement among the SQ grades (0.26–0.45). Moreover, the accuracy assessment values varied from 39–49% for TDS and MDS. Previous studies [11,16] showed that the RF model accurately predicted SQ grades in semiarid and humid regions.

3.5. Prediction Map of SQ Grades

The different SQIs maps were prepared using the RF model based on the DSM framework (Figure 3). As previously discussed, all SQIs are classified into five classes (Table 6). According to the predicted maps, the dominant class of SQ is moderate (III) and high (II) exception SQI_n-L-MDS showed the dominant SQ grade of IV (low) (Figure 3A). Only a small area is located on grades I (very high) and V (very low) SQ on predicted maps based on the L scoring system (Figure 3A). Representing maps of all SQIs showed the lowest proportions of very high and high SQ classes in the north of the area where pistachio and barley land covers were located, while pomegranate and saffron have the lowest proportions of very low SQ (grade V) (Table 7 and Figure 3). This is because higher soil fertility parameters (i.e., SOC, TN, P_av, K_av, and SP) and lower salinity and alkalinity correspond to better SQ availability through the pomegranate and saffron land uses. Previous studies [7,16,49,53,54] have reported that areas with optimal soil properties are known as high and very high SQ areas. The EC, Na_aq, and SAR values point out the crucial effects on crop yields, and also salinity and high Na_aq contents can have an unfavorable impact on soil aeration and physical properties, including soil structure and infiltration [49], as a consequence of land degradation.

On the predicted maps based on the NL scoring system (Figure 3B) in SQI_w and SQI_n, the areas are commonly categorized as grade I, II, and V. It should be noted that the very low grades were located in the north of the area with the cultivation of pistachio and barley. Thus, these findings confirmed the effects of landform (i.e., playa margins) environmental covariates and farmer management on the SQ of the studied region.

The results of the evaluation of the RF model for predicted maps are represented in

Table 8. Validation criteria of SQIs maps and random forest (RF) model result for SQI predictions.

SQI	RMSE	R2	MAE	RMSE + SD	R2 + SD	MAE + SD	mtry
SQIw-L-TDS	0.049	0.373	0.038	0.005	0.099	0.004	24
SQIw-NL-TDS	0.111	0.373	0.092	0.013	0.131	0.011	54
SQIw-L-MDS	0.075	0.174	0.061	0.009	0.114	0.008	31
SQIw-NL-MDS	0.146	0.273	0.120	0.012	0.116	0.009	39
SQIn-L-TDS	0.033	0.391	0.026	0.003	0.114	0.002	31
SQIn-NL-TDS	0.074	0.393	0.062	0.008	0.111	0.007	31
SQIn-NL-MDS	0.088	0.254	0.072	0.007	0.102	0.006	9
SQIn-L-MDS	0.051	0.152	0.040	0.004	0.073	0.003	9

SQI_w—weighted additive soil quality index; SQI_n—nemoro soil quality index; TDS—total data set; MDS—minimum data set; L—linear and NL—nonlinear.

. Predicted maps of SQI_n-L-TDS and SQI_n-NL-TDS had the best performance with R² = 0.39. In contrast, the lowest predicted model (R² = 0.15) was found for SQI_n-L-MDS (Table 8). The important point is that predicted map SQI_n-L-MDS is different from other maps and mostly predicted as a low (IV) class due to the RF model’s low accuracy. Many studies have used the RF model to predict soil properties and reported high prediction performances of RF [2,16,21,23]. In addition, the results of this research demonstrated the high capability of the RF model for the SQI prediction map. Geomorphological conditions of the study area have affected many soil properties [6] and consequently SQI in the region. The RF model had a relatively accurate estimate of the SQI. Therefore, it is suggested to use machine learning techniques such as RF and auxiliary data such as geomorphological maps, topographic attributes, and satellite images to map SQI.

3.6. Environmental Variable Importance

As shown in Figure 7, the most important environmental covariates were obtained from the relative variable importance analysis of the RF model. The results showed that the covariates, including wind effect, multiresolution valley bottom flatness index (MRVBF), valley depth, and groundwater quality parameters (i.e., pH, SAR, HCO₃⁻), were considered as the most important factors to predict the spatial SQIs variability based on TDS and MDS (Figure 7). Moreover, the results indicated that in the SQI_{n_L_TDS} map, groundwater quality parameters were the most important and made a higher contribution in predicting SQI maps. However, in the SQI_{n_NL_MDS} map, the topographic attributes were the most important covariate and showed a higher contribution. In line with our results, the previous studies showed that topography had an important role in the spatial variation of soil properties and SQ in different regions of Iran [2,11,16,23]. These researchers reported the role of topography attributes in soil erosional and depositional processes of water and particles [55,56,57]. Moreover, the distribution of soil salinity and alkalinity is mostly controlled by topography attributes. The lower grades of SQ were identified in the lowest elevation (playa margins or northern of the study area, Figure 1) with the highest soil salinity and alkalinity levels. Of particular importance is the wind effect in the study area, which should be noted due to the high salinity and lack of fertile lands as half of the study area lacks vegetation, meaning the area is prone to wind erosion.

On the other hand, this study identified the effects of irrigation water quality parameters in SQI map predictions and consequently soil salinization and provided a scientific basis for improving the Bajestan ecosystem maintaining the healthy development of agriculture. Several studies have assessed the driving factors of groundwater quality on soil salinization, e.g., [2,10,58]. They reported that frequent use of inadequate irrigation water could intensify soil salt aggregation in agricultural lands that increased soil salinity and decreased SQ and shallow groundwater, seriously threatening the long-term sustainability of irrigated agriculture. Thus, the Bajestan area requires considerable attention, especially in the north area and pistachio orchards, as groundwater quality is reduced from south to north [6]. Furthermore, saline groundwater used for irrigation in the cultivated lands has increased salt concentration in the root zone, and because of low rainfall in the area, there is insufficient leaching, and long-term use of low-quality water in the area can lead to increased salt concentration in soil and decreased SQ.

4. Conclusions

In this study, two SQIs with two scoring systems were used to assess SQ in an agroecosystem. The results of this study illustrated that the lands under saffron and pomegranate cultivations have higher SQ than pistachio and barley, as clearly showed the SQIs in TDS and MDS. According to the results obtained between the two data sets for all SQIs, the MDS performed as a reliable, fast, and economically appropriate solution to select the minimum effective soil properties in the study area can be useful. The application of this method reduces the reproducibility effect of similarly correlated properties and can display the information in other parameters as a selected set. The RF model showed promising results in the distribution of predicted SQ grades.

The SQ’s spatial distribution maps revealed that the study area’s SQ was mainly moderate (III) and high (II). Moreover, the SQI_n could provide a better estimation than SQI_w in the L and NL scoring system of TDS. The results showed that the topographic attributes and groundwater quality parameters had the highest effect on SQI distribution in the study area. Therefore, to plan for sustainable agriculture and cleaner production, appropriate management and protection measures should be taken in the lands under pistachio and barley cultivation due to the low SQ situation in the region’s northern areas. Furthermore, proper irrigation is necessary for arid areas to leach the salts from the root zone and increase water efficiency. The results proved that SQ assessment has good efficiency for evaluating soil health status in different locations of areas for decision making on cultivation patterns. Consequently, any management and crop that reduces the soil surface vegetation and uses inadequate irrigation water could reduce the SQ and suitability. Overall, sustainable soil management in the Bajestan area requires the development of efficient management guidelines and the provision of appropriate management methods in the field of fertilization, tillage, and land irrigation to farmers.