Assessment of Agricultural Soil Quality Indices Using Mechanistic Models ()
1. Introduction
The global security of soil and water resources has been facing numerous challenges and concerns associated with the rapid increase in the human population, inefficient use of soil and water resources, intensive agricultural activities, and changes in land use [1] - [6] . All these human activities have led to the overexploitation of the ecosystem [7] [8] , threatening soil quality, resulting in soil degradation and a reduction in its capacity to sustain plant and animal production within natural and managed ecosystems. Consequently, soil degradation has raised concerns regarding biogeochemical cycles, climate change, and ecosystem sustainability [9] [10] [11] . Therefore, maintaining soil quality is critical to sustain and strengthen ecosystem services, ensure the high nutritional value of the food chain [12] [13] and promote the growth and productivity of agricultural products [14] . Furthermore, poor agricultural practices ultimately negatively affect human health. The accumulation of heavy metals, such as zinc, has been intensively studied for its adverse impacts, including nausea, abdominal pain, anemia, dizziness, lethargy, and gastrointestinal effects [15] [16] [17] .
Soil management strategies play a crucial role in governing soil quality [18] . The reduction in soil quality arises from cultivated lands shrinking due to improper land use and overexploitation of landscapes [19] [20] [21] . Surface soil erosion stands out as one of the main factors contributing to the reduction in soil quality [22] . Soil erosion primarily leads to a decrease in soil organic matter, eutrophication, and nonpoint source pollution [23] [24] [25] . Burning crop residues also contribute to soil quality reduction as it has the potential to decrease microbial biomass and increase nutrient immobilization [24] [25] . It has been reported adverse impacts of deforestation and agricultural practices on soil quality [26] . Earlier studies [27] [28] employed fuzzy methods and spatially multivariate techniques to estimate desertification and degradation rates and evaluate soil quality in Turkey.
Research on soil quality demands more attention to the soil’s biological, chemical, and physical properties simultaneously [9] . However, soil quality cannot be directly calculated and must be inferred from standard soil quality indicators. Soil quality indicators are measurable soil characteristics that contribute to soil capacity for crop production or environmental performance and are sensitive to land use change, management, or agronomy practices [9] . Several tools or packages have been identified for assessing soil quality, including soil card design and farm packages, the soil quality index [29] , and geostatistical methods [30] . Among these packages, the soil quality index is the most widely used method [29] because it is easy to use and particularly useful for incorporating temporal and spatial information [31] . Soil quality is dynamic and depends on the interactions among soil factors, making it challenging to accurately identify soil factors for management decisions. Among the methods for soil factor selection, the total data set (TDS) and minimum data set (MDS) have mainly been used in soil quality assessment studies [32] [33] . In TDS, a set of total soil factors is measured based on the analyzed soil characteristics, while MDS only relies on some selected soil factors based on the correlation matrix between the indicators, which decreases the complexity of the data and also enhances pertinent information through correlation of the soil indices [34] . Determining the soil quality index as an indirect method is based on integrating dimensionless factors (extracted by normalization) with their weights. Consistent with that, geostatistical methods can further provide efficient and reliable tools to increase the number of measurement points in unsampled locations [35] . Among the geostatistical methods, kriging interpolation is considered ideal and widely used among other interpolation methods [36] . This method estimates the values of indeterminate points using the values of specific points. The kriging technique results are more effective when the variables have a normal distribution [35] [37] . As a technical method, a geographic information system (GIS) combines geographical features using a data table to analyze, map, and evaluate problems. The main strength of GIS is the use of spatial and statistical methods to analyze the distribution and geography of information.
Many farms in the north of Iran have been abandoned due to high salinity and the lack of a suitable reference for making fertilizer recommendations. Studies conducted so far have only addressed the physicochemical and morphological properties of the Golestan province, without conducting zone-to-zone soil physicochemical evaluation [38] . Furthermore, a simple and suitable method for assessing the soil quality of agricultural areas has not been provided in Golestan province to assist farmers in practical management. Given the lack of appropriate studies on soil quality and management decisions, GIS-based qualitative research presents an opportunity to help growers and policymakers make more reliable decisions to reduce soil erosion and mismanagement practices.
This study aims to assess agricultural soil quality in the northern region of Golestan province, Iran. It employs two soil quality assessment indicators, IQI and NQI, along with two datasets, TDS and MDS. Additionally, the study utilizes GIS for the spatial display of these indicators. The objective is to identify and select the most appropriate soil quality indicators for this region using statistical methods, linear relationships, and match analysis. Furthermore, the findings of this study could serve as an initial approach for future studies in developing nations like Iran, where the assessment of indicators needs to be cost-effective and feasible with minimal infrastructure to ensure widespread adoption.
2. Materials & Methods
2.1. Description of the Experimental Site
The study area is located in the eastern part of Kolijeh village in Golestan Province, Iran. Kolijeh village is situated in the DashliBorun district, at 55˚40' east longitude and 37˚57' north latitude (see Figure 1). The altitude from the sea level is approximately 70 meters, with an annual rainfall of 170 mm and an average surface evaporation of 1900 mm. The mean annual temperature, along with the
![]()
Figure 1. Study area and sampling point.
average minimum and maximum temperatures, are recorded as 16.8˚C, −1˚C, and 36.6˚C, respectively. The majority of rainfall occurs during the spring season, typically from late March to late May. The slope of the study area is approximately 2%, and the general slope of the site is towards the north. The soil in this area is classified as inceptisol according to the USDA soil classification system. Agricultural activities in the area are predominantly focused on cultivating wheat, barley, and blue cotton crops.
2.2. Soil Sampling, Management, and Measurement of Characteristics
Soil samples were collected twice in 2022. The first sampling took place in early fall, after harvesting the fall wheat, while the second sampling occurred just before the start of the new growing season in spring. As shown in Figure 1, the sampling network distribution in the study area was designed with 750 × 750 m2 networks, with the center of each serving as the sampling point. A total of 89 soil samples were collected, with the soil sampling depth ranging from 0 to 20 cm. The geographical coordinates of these samples were recorded using a handheld GPS device. The collected soil samples were subsequently laboratory-dried and passed through a 2 mm sieve for physical and chemical analysis. Various soil properties were measured, collectively referred to as Total Data Set (TDS), including pH, electrical conductivity (EC), organic carbon (OC), percentage of calcium carbonate equivalent (CCE), sand, silt, clay, nitrogen (N), potassium (K), available phosphorus (P), sulfur (S), iron (Fe), zinc (Zn), manganese (Mn), and copper (Cu) (as listed in Table 1). Essential nutrients such as N, P, Zn, and Fe are crucial for ensuring optimal plant growth, maturation, and overall yield. Any deficiency in these nutrients can lead to abnormalities in plant development and a reduction in productivity. Therefore, our focus was primarily on these elements to underscore their significance in the soil-crop-environment nexus.
3. Evaluation of Soil Quality Index
3.1. Scoring Characteristics
Standard Scoring Functions (SSF, see Table 2) were employed to assess the characteristics [46] , with the values of each index being converted to scales 0 and 1. Three functions were utilized based on the sensitivity of the soil quality characteristics:
1) “More is better”: This function was applied to score soil organic carbon characteristics, high consumption, and low consumption elements that play a crucial role in soil fertility, soil water retention, and accessibility.
2) “Less is better”: This function was used to evaluate soil acidity and lime properties because high concentrations of these properties in arid and semiarid regions can reduce plant access to nutrients [39] .
3) “Appropriate limit”: This function was employed for determining the EC level. Considering the optimal pH value of 7 and the soil pH variability over 7 in all measured points, pH was categorized under the “Less is better” group instead of the optimal category. Soil EC, which ranged between 1.01 and 10 dS∙m−1, was placed in the “Appropriate limit” group (0.2 - 2). For points less than 0.2, the “More is better” function was applied, while for points higher than 2, the “Less is better” function was utilized. Additionally, since there is no information available on the upper and lower thresholds in the study area, the minimum and maximum values of the variables were considered the low and high thresholds, respectively. Each characteristic’s F(X) score was calculated using these values in the standard scoring functions.
![]()
Table 1. Soil physical and chemical analyses for selected indicators.
![]()
Table 2. Standard scoring functions and measured characteristics in Golestan province.
FTa means function type; LB means less is better; MB means more is better; AL means appropriate limit. SSFb stands for Standard Scoring Function. X means the characteristic value, F(X) means the value of the characteristics between 0.1 and 1, and L and U are the low and upper threshold values, respectively.
3.2. Selection of the Minimum Data Set
MDS was derived through principal component analysis (PCA) using MINITAB software version 17.0 to identify the most influential parameters from TDS and to minimize measurement costs. In this method, a correlation matrix between the variables was constructed, followed by the calculation of specific vectors and values within this matrix. The majority of components (PCs) with eigenvalues greater than or equal to one (Eigenvalues ≥ 1) were utilized for selecting MDS [47] . However, in certain cases, PCs with eigenvalues less than one were also considered to enhance the variance explained by the model. Once the required number of PCs was determined, variables with the highest weights (specific vectors) in each PC were identified. Initially, the variable with the highest weight in the desired PC was selected, and other variables with weights exceeding 90% of that weight were designated as the most critical variables of that PC. It is crucial for the selected variables to exhibit significant correlations; the variable with the highest weight represents the group. However, if there is no correlation between a variable and others, more than one variable may be chosen in each PC as MDS. The plots were generated using the ggplot2 package in RStudio version 2023.06.0 [48] , and the map visualization was conducted using ArcGIS 10.3.
3.3. Determination of Weight Coefficient and Soil Quality Indicators
Principal factor analysis was employed to determine the weighting coefficient for each variable. Through this analysis, the common variance of the variables for both TDS and MDS was extracted. The weight coefficient for each variable was calculated from the ratio of the common variance of each variable to the total common variance of all variables, using SPSS version 16. Using the following equation, the value of the IQI index is determined at each point (1).
(1)
where: Wi is each indicator’s weight, Ni is each attribute’s score, and i is the number of variables or attributes.
And the Nemoro Soil Quality Index (NQI) Equation (2).
(2)
Pave is the average score of the indicators (derived from the standard scoring function) at each point, Pmin is the minimum score of the selected indicators at each point, and n is the number of indicators.
3.4. Soil Quality Classes
Jenks’s optimization method was employed to group the values of each index into 5 classes [49] . Jenks’s method operates to minimize the variance within the classes while maximizing the variance between the classes. Soil quality characteristics exhibit dynamic, changeable, and spatially varying attributes. Understanding these variations is crucial for implementing site-specific management practices in precision agriculture. Therefore, one parameter expressing these characteristics can be considered a variable, and its variability can be studied through statistical methods. Recent advances in statistical theory have facilitated the quantification of spatial relationships between samples at smaller distances. The foundation for this advancement lies in the theory of regional variables known as Geostatistics. Geostatistics, a branch of applied statistical science, offers a wide range of statistical estimators to predict the desired property at non-sampled locations using information from sampled points [30] . In this research, the kriging method was employed to analyze the structure of soil quality changes, predict trends, and estimate properties.
3.5. Index Evaluation
After determining the most appropriate interpolation method to evaluate the agreement of the indicators and their estimation accuracy in MDS compared to TDS, conformity analysis was conducted using direct comparison, kappa coefficient, correlation analysis, and regression. Direct comparison assesses the number of points with the same soil quality class. The kappa index measures the inter-rater or inter-commentator qualitative index, which outperforms simple agreement calculations. In other words, it provides an index in which the effect of chance agreement is adjusted, thus always ranging between −1 and +1. Values close to plus one indicates high agreement, while those close to minus one indicate disagreement.
4. Results
4.1. Soil Quality Assessment Based on TDS
Table 3 presents the standard deviation, mean value, coefficient of variation, minimum, maximum, range, skewness, and kurtosis of the 14 indicators measured at each sampling point. The highest levels of K, Mn, and Cu were obtained from the zones under investigation recorded in the study area, while the lowest levels of N, P, Zn, and Fe were observed in the area. Significant changes in soil pH and CCE content were recorded throughout the study area, while low electrical conductivity levels were observed only in soils near the Atrak River at 7 Ds m-1.Pearson’s correlation coefficient (Table 4) for all indicators showed that the most potent synergistic correlation exists between N and OC (r > 0.80), followed by Fe and Cu (r > 0.68). Conversely, the most robust antagonistic relationship is between clay and silt (r < −0.75), followed by silt and Cu (r < −0.46). Based on the estimated communality analysis, the values of Zn, Fe, and P received the highest weight (average weight is 0.08); in contrast, soil reaction parameters and calcium carbonate equivalent received the lowest weight (average weight is equal to 0.04). According to Table 5, IQI and NQI indices for the TDS were grouped into five classes and had a range of values between <0.678 (very high) to <0.545 (very low) for the IQI index and >0.472 (very high) to <0.370 (very low) for the NQI. Five model performances (Circular, spherical, tetrahedron, exponential, and Gaussian) (Table 6) were compared to fit the experimental half-view. The root mean square error (RMSE) was used to compare the fitted models to the experimental half-shift. Based on the mutual evaluation, the model fitted to the half of the experimental view change was an exponential model with the root of the standardized square mean close to one. The maps derived from the IQI model (Figure 2) revealed that the soil quality of the study area was predominantly very low (grade 5), covering an area of 3096 km2, which accounted for 82% of the study site. A small portion of the study area exhibited low-quality grade (4) and medium grade (3), constituting 11% (473 km2) and 5% (215 km2) of the total area, respectively. Similarly, the soil quality map obtained from the NQI model yielded comparable results to the IQI model, with the entire study area demonstrating very low-quality grade (grade 5) (Figure 2). In both instances, the soil quality exhibited a decline from south to north (Figure 2).
![]()
Table 3. The mean, standard deviation, coefficient of variation, minimum, maximum, and range of indicators.
EC: electrical conductivity.
![]()
Table 4. Pearson’s correlation matrix between soil properties.
![]()
Table 5. Classification of soil quality levels for the studied indicators.
*I (Very High), II (High), III (Moderate), IV (Low), V (Very Low).
![]()
Table 6. Experimental models based on experimental variogram for each of the soil quality assessment indicators.
a: Root Mean Square Error.
![]()
Figure 2. Distribution of soil quality levels with different indicators in Golestan province.
4.2. Soil Quality Assessment Based on MDS
The most effective TDS parameters were determined by extracting the MDS using principal component analysis (PCA). Five factors or PCs were extracted out of 14 variables, with the range of eigenvalues of these components (PCs) being between 4.3 to 1 (Table 7). According to the analyses, five variables, namely EC, P, Fe, silt, and clay, were selected as MDS, collectively explaining 76% of the variance changes. Based on the results of the estimated commonality analysis of variance, the lowest value (0.122 value of common variance) was observed for the EC of the soil, which received the lowest weight (0.038 value weight). Conversely, the variables Fe, P, clay, and silt exhibited a high contribution in defining soil quality (Table 7).
Similar to the TDS method, in the MDS approach, the soil quality of the study area was categorized into five classes, with a grade ranging between 0.675 - 0.438 for the IQI index and a range of 0.427 - 0.264 for the NQI index. Through the cross-evaluation method, the exponential model fitted to the semi-experimental exponential change was selected as the optimal model for both IQI and NQI indices (Table 5). The soil quality map based on the NQI index indicated that 93% of the study area exhibited a very low-quality level (4004 km2), while only 7% of the study area (308 km2) demonstrated an average soil quality grade. In contrast, according to the IQI index, 69% of the study area was classified as having very low and low quality (2975 km2), while another 17% of the study area was identified as possessing medium (603 km2) and high (733 km2) soil quality, respectively. Similar to TDS, maps generated from the two indicators applied to MDS illustrated an increase in soil quality from south to north of agricultural land.
![]()
Table 7. Results of principal component analysis (PCA) of soil quality indicators of agricultural lands in Golestan province.
a: Main components; b: Selected factors in MDS.
4.3. Validation of Models
Direct analysis for 89 sampling points between TDS and MDS-based soil quality indices revealed a 95% and 90% agreement, respectively. The agreement between TDS and MDS characteristics for NQI and IQI models was 32% and 43%, respectively. Similar results were obtained by applying kappa-specific analysis. The kappa value between IQI and NQI using the TDS method was 0.6, indicating significant agreement, while in the MDS method, a value of 0.1 was obtained, indicating weak agreement. For IQI and NQI models, the kappa coefficient values were 0.15 and 0.111 (P < 0.05) between TDS and MDS, respectively, indicating low agreement. Linear relationships between different methods showed a higher correlation coefficient (R2 = 0.43) when the IQI model was used (Figure 3).
![]()
Figure 3. Linear regression between TDS and MDS methods for NQI and IQI models.
5. Discussion
5.1. Soil Quality and Characteristics
As Table 8 represents, all soil samples contained a certain amount of organic carbon (<5%), which is typical for arid and semiarid soils [50] [51] . The soil EC values in the study area were classified into three categories: non-saline soils, low salinity, and soils near the starting point of the Atrak River, where salinity levels were close to seven. High levels of equivalent calcium carbonate (CCE > 15%) and soil pH (>7) in agricultural lands indicate alkaline soil types. Similar results were obtained in previous studies of arable soils in Golestan province.
Due to extensive agricultural operations, low amounts of trace elements (Fe and Zn), high consumption of N and P, and soil organic carbon have contributed to the decline in agricultural land quality. Factors such as failure to use fertilizer based on soil tests, burning of plant residues after harvest, and lack of proper drainage systems can be attributed to the main reasons for reducing soil quality in the region.
On the other hand, in calcareous soils of Iran, despite the abundance of some nutrients such as P, Zn, Fe, and Cu, the high pH and abundance of calcium result in the available forms of these elements being less than necessary for plant growth. Therefore, plant productivity is generally limited by nutrient deficiencies [52] [53] .
Soil quality indicators (IQITDS, IQIMDS, NQITDS, NQIMDS) indicated that the majority of agricultural lands in the study area (approximately 71%) exhibited very low soil quality (grade 4 and grade 5), with only about 5% of the lands
![]()
Table 8. Results of estimating common values of variance and weight of each indicator for both TDS and MDS methods.
*Communality.
showing very high soil quality (grade 1). In the TDS method, IQI and NQI indices revealed low and very low-quality soils for most agricultural lands (Grade 4 and Grade 5, Figure 2), despite the observation of the highest EC values (up to 7.4) in this region. This phenomenon could be attributed to the proximity of these areas to the Atrak River and fluctuations in groundwater levels. The use of water extracted from the groundwater contributes to salt accumulation and increases the soil surface’s EC. Soil salinity and its adverse effects on agriculture, coupled with uncontrolled agricultural practices, may lead to a reduction in soil quality.
A high percentage of calcium carbonate (>15%) was detected in the study area. This level of calcium carbonate, in conjunction with high pH (>7.5), restricts the uptake of certain elements such as Fe, Zn, and P by plants [54] . In the MDS method, the maps derived from IQI and NQI indices depicted a decline in soil quality from south to north (Figure 2), correlating with the amount of clay and subsequently, a decrease in soil organic carbon. In arid and semiarid regions, the quantity of organic matter is influenced by soil properties such as texture, clay content, minerals, EC, and soil acidity [55] . Soil texture serves as a fundamental qualitative and physical property [56] . Numerous studies have noted that soil texture significantly impacts the properties of soil organic matter, with soils possessing high clay content exhibiting the highest levels of organic matter and the most effective physical and chemical protection of soil organic matter [55] .
All soil samples had an average percentage of clay and silt of 20% and 38%, respectively (Table 3). The low level of fine soil particles, such as clay, in the study area has hindered the maintenance and stabilization of organic carbon. Additionally, high EC adversely affects the soil microbiome, diminishes soil quality, and ultimately constraints crop production [24] .
5.2. Indicators and Characteristic Methods
In various studies, PCA has been explored as a tool for data reduction [31] . This method not only selects characteristics that better represent soil quality but also saves time and reduces the costs associated with laboratory analysis. Among the characteristics selected in the MDS method, previous studies have identified EC, soil P, percentage of clay particles, and soil Fe as crucial indicators [56] [57] .
In statistics, the Kriging estimator holds significance as one of the most crucial linear estimators due to its minimal estimation variance and lack of systematic error. The zero mean estimation error is a primary requirement for Kriging [57] . Kriging has been applied in various spatial analysis topics, including the assessment of soil quality in agricultural lands, as demonstrated in this study.
Conformity analysis between indicators and characteristic methods yielded low values, primarily due to the absence of a standardized soil classification method and critical values specific to the study area, particularly in Iran. These findings align with previous studies conducted in Iran’s agricultural lands [56] .
In our study, analysis using kappa agreement and direct comparison revealed higher IQI index values compared to the NQI index for both TDS and MDS methods. Regression and correlation coefficients further indicated the superior estimation of soil quality through the application of the IQI index over the NQI index. This difference may be attributed to the characteristic weights, which serve as a distinct factor in the IQI model, as opposed to the NQI model, which prioritizes minimal crop production.
Although kappa analysis demonstrated a higher value for TDS compared to MDS, our results underscored a consistent trend in soil quality assessment with the application of the IQIMDS model. These findings were consistent with those of [31] . in their study of soil characteristics in agricultural lands in Qazvin province, Iran.
6. Conclusions
This study compares two sets of characteristics, TDS and MDS, along with two different indicators, IQI and NQI, to assess the soil quality of agricultural lands in the north of Golestan province, Iran. All measured characteristics were deemed suitable for soil quality assessment. However, among these characteristics, only EC, P, Fe, clay, and silt had the highest impact on soil quality and were included in the MDS through PCA.
Direct analysis of 89 sampling points between TDS and MDS-based soil quality indices revealed a 95% and 90% agreement, respectively. The associated R2 for TDS and MDS characteristics for NQI and IQI models were 32% and 43%, respectively. These values indicate that most soils in the studied region exhibit low quality across all studied indicators (IQITDS, IQIMDS, NQITDS, NQIMDS). While the match analysis showed good agreement for IQITDS, a favorable outcome was obtained for IQIMDS. Consequently, the IQIMDS method can be considered a suitable tool for developing a quantitative soil quality assessment method.
Future research should focus on further delineating lands into specific zones and applying tailored management systems according to the unique needs of each zone. These results can serve as a starting point for future researchers to develop improved soil assessment tools and explore other soil quality models. Modeling soil characteristics, albeit complex, holds promise for effective farm management, predicting agricultural productivity, and closing the yield gap.
Authors Contribution
Nikpey M.*: Conceptualization, Investigation, Methodology, Data Curation, Data Visualization, Validation, Writing-original Draft, Robatjazi J.*: Conceptualization, Investigation, Methodology, Data curation, Writing-original Draft, Writing-review, Data Visualization, Validation, Supervision of The Research, *These authors contributed equally, Garmehei Sh.: Investigation, Data Curation, Writing-review, Hendra Gonsalve W Lasar: Writing-reviewing, editing, Nguyen Khoi Nghia: Writing-reviewing, editing, Kwadwo Agyei B.: Review, editing. All the authors contributed to the editing of the final draft.
Acknowledgments
The authors would thank the reviewers for their time and constructive feedback in this paper.
Funding Source
This study received no funding.