# Evaluation of Soil Quality in Arid Western Fringes of the Nile Delta for Sustainable Agriculture

Authors: Ahmed M. Saleh; Mohamed M. Elsharkawy; Mohamed A. E. AbdelRahman; Sayed M. Arafat
Journal: Applied and Environmental Soil Science
 (2021)

---

## Abstract

Egypt is currently witnessing an extensive desert greening plan with a target of adding one and a half million feddans to the agricultural area. The present study evaluates the soil quality in the western desert fringes of the Nile Delta using three indicator datasets, which involve the total dataset (TDS), the minimum dataset (MDS), and the expert dataset (EDS). Three quality index models are included: the Additive Soil Quality Index (SQI-A), the Weighted Additive Soil Quality Index (SQI-W), and the Nemoro Soil Quality Index (SQI-N). Linear and nonlinear scoring functions are evaluated for scoring soil and terrain indicators. Thirteen soil quality indicators and three terrain indicators were measured in 397 sampling sites for soil quality evaluation. Factor analyses determined five soil and terrain indicators for the minimum dataset and their associated weights. The linear scoring functions reflected the soil system functions more than nonlinear scoring functions. Soil quality estimation by the minimum dataset (MDS) and Weighted Additive Soil Quality Index (SQI-W) is more sensitive than that by SQI-A and SQI-N quality models to explain soil quality indicators. The moderate soil quality grade is the largest quality grade in the studied area. The minimum dataset of soil quality indicators could assist in reducing time and cost of evaluating soil quality and monitoring the temporal changes in soil quality of the region due to the increased agricultural development.


---

## Body

## 1. Introduction

Soil is considered one of the most important natural resources for countries where the foundations of human life are mainly based on their suitability for agriculture as well as other anthropogenic activities [1]. In the Eastern Mediterranean, the Egyptian civilization has sustained in the fertile Nile Valley and the Nile Delta, in which the elements of sustainability are available [2]. During the past century, there has been a serious deterioration in the productive capacity of more than 10% of the world’s lands, which increased the interest in assessing soil health at the regional level [3]. Assessment of soil quality is necessary to assist farmers in evaluating the effects of their management decisions on soil productivity [4, 5]. The study in [6] emphasized that soil qualitative evaluation, despite its simplicity and speed, requires very high experience. Accurate assessment of soil quality is complex due to the lack of consistency in soil characteristics [7–10] and the different soil management practices such as adding organic and mineral fertilizers, herbicides, and pesticides, which make the soil assessment more difficult [9]. Soil quality indicators are evaluated using different descriptive and quantitative approaches [6, 11–13]. The concept of soil quality is more comprehensive than the reductive approach of measuring a single indicator [14]. Basically, soil quality indices have three component targets as the environmental quality, agronomic sustainability, and socioeconomic feasibility [15]. Therefore, the estimation of soil quality indices is a complicated process and quite a hard mission [16]. Recently, several techniques have been successfully utilized for evaluating soil quality such as spectroscopy fingerprint [17], the Soil Management Assessment Framework (SMAF) [18], and spatial variability of physicochemical soil properties [19]. Most researchers estimated soil quality using only one technique with some exclusion [17–21]. Principal component analysis (PCA) and linear regression coefficients are also widely used for minimizing the number of soil variables included in estimating soil quality [21]. Recent studies attempted to develop appropriate algorithms and scoring functions for determining soil quality [22] by comparing multiple soil quality scoring techniques [15, 23, 24].To overcome the financial constraints, the amount of soil quality indicators needs to be reduced to a minimum dataset (MDS) [25, 26]. The concept of minimum datasets was proposed in [27], which is the lowest set of soil quality indicators required to measure or evaluate soil quality to demonstrate the soil’s capacity for sustainable agricultural production. Recent soil quality studies preferred the minimum dataset (MDS) approach [28–30], where the MDS was favorable for detecting yield variations than the total dataset approach [31]. Additionally, the MDS approach is an effective tool for identifying the status of soil nutrients at the regional level [32]. The limitation factors of soil quality vary based on the soil topography, land use, climatic zone, the regional soil ecosystem, and the soil geological units [33, 34].By establishing the MDS, it is easy to select the most effective indicators to determine soil quality determined based on its ability to estimate soil productivity and soil stability, and this indicator has been widely used [35–37]. Based on MDS indicators, the abundance of soil data can be minimized [10, 35]. Furthermore, the weight of chosen factors can be estimated while establishing the minimum datasets, thus minimizing the internal influence [37]. The MDS can be determined by linear or multiple regression factor analysis [38], multivariate statistical techniques [9, 37], discriminant analysis, and score function [34]. Since factor analysis can minimize abundant input data in the original soil datasets, it is widely utilized in the limitation of the MDS [28, 39]. Factor analysis statistical techniques can identify the most important indicators discriminating soil quality in combined tillage, fertilization, and crop rotation treatments [28]. In Egypt, recent researchers evaluated soil quality in the Nile Delta and El-Fayoum depression [40–42].As the calculation of soil quality is hard [4], an urgent need exists for improving the simple and reliable quantitative evaluation of soil quality through comparison of various available assessment strategies and methods. Thus, this study attempts to examine the soil quality in the arid region along the desert western fringes of the Nile Delta by (i) assessing the soil quality using three types of indicator datasets (TDS, MDS, and EDS), three types of linear scoring functions, two types of nonlinear scoring functions, and three soil quality index models (additive, weighted additive, and Nemoro quality indices); (ii) suggesting the most suitable indicator method and soil quality index model for the studied region using sensitivity analysis and linear relationships. The findings of the study have a significant value for assessing the soil quality in the rapidly developing area for sustainable management.


## 2. Materials and Methods

### 2.1. Field Description of the Study Area

The evaluated area is located in Wadi El Natrun district, El Beheira Governorate, at the western fringes of the Nile Delta between 29°54′00″E–30°20′00″E and latitudes of 30°22′00″N–30°00′00″N (Figure1). The evaluated area covers an area of 1600 km2 (380000 feddan), with a landscape consisting of Wadi terraces and Wadi depressions with longitudinal sand dunes at the southern edge of the area [43]. The annual evaporation rate is about 114.3 mm/y, the average air temperature is 21°C, and the area has a rarely annual precipitation rate with an average of 41.4 mm/y [44] with an erratic distribution. Based on the work in [45], the soil temperature regime is “Thermic” and the moisture regime is “Torric.” Two soil orders are dominant in the area, Entisols and Aridisols. The evaluated area has been reclaimed during the last two decades and is cultivated with diverse kinds of fruits and vegetables [46]. This kind of cultivation is usually associated with an increase of underground water irrigation, and the ongoing climate change put a threat on the soil quality and sustainable agriculture that reduces the economic value of the lands [47, 48]. Additionally, the excessive extraction of groundwater for irrigation resulted in land subsidence, which could also affect the integrity of the soil profiles in the area [49].Figure 1
Location of the study area and soil sampling sites.


### 2.2. Data Collection

#### 2.2.1. Indicator Sampling Design

The required samples size needed for soil quality assessment was determined by using the binomial probability algorithm [50]. 397 sampling sites were chosen. A Proportionate Stratified Random (PSR) sampling technique was used as the sampling type for soil and terrain sampling based on land cover and land use distribution in the evaluated area. The sampling locations were georeferenced (Figure 1) using the Global Positioning System (GPS).


#### 2.2.2. Sampling and Analyses

The soil sampling was carried out according to the sampling design. However, some sampling sites were displaced by other sampling sites near the designed sampling locations due to accessibility difficulties, and some certain areas were currently restricted for access, also the presence of elongated sand dunes at the south of the study area. These accessibility limitations resulted in ununiform spatial sampling distribution.Disturbed and undisturbed samples were collected to the depth of the root zone (30 cm) and described according to the work in [51]. Soil samples were crumbled gently by hand without root material and were air dried, grounded, and sieved through a 2 mm sieve. Soil chemical analysis, including pH in 1 : 2.5 soil-water suspension, electrical conductivity (EC) in soil paste extract, cation exchange capacity (CEC), exchangeable sodium percentage (ESP), organic matter (OM), and calcium carbonate were determined according to standard methods of [52]. Soil physical analyses, including particle size distribution using the pipette method, bulk density (BD) using a core method, available water holding capacity (AWHC), and field capacity (FC), were performed according to the work in [53]. Soil packing density (PD) was estimated as a composite index of soil bulk density and soil clay content according to the work in [54].Terrain analyses were performed on a shuttle radar topography mission (SRTM) 1-arc-second v.30 image [55] to obtain some surface parametric information (elevation, slope, relief intensity, and topographic witness index) [56]. The DEM image was firstly sink filled using depression filling algorithm [57]. The slope (Slp) and topographic witness index (TWI) were extracted [58, 59]. The extraction of relief intensity (Rlf) was based on the algorithm described in [60].The Normalized Difference Vegetation Index (NDVI) and elevation data were used as exterior environmental variables in the process of representative indicator selection and redundancy reduction. NDVI data were retrieved from the Landsat Data Continuity Mission (LDCM) sensor multiband dataset using the NDVI equation ([61, 62]).All of the datasets were projected into the WGS84-based Universal Transverse Mercator (UTM) orthographic projection coordinate system (EPSG 32636) and resampled to a 30 m spatial resolution.


### 2.3. Indicator Selection and Scaling

Three datasets of quality indicators were chosen, the total dataset (TDS), minimum dataset (MDS), and expert opinion dataset (EDS). The total dataset (TDS) includes all soil and terrain indicators for soil quality index development. The soil and terrain indicators for the minimum dataset were determined according to multivariate factor analysis. Factor analysis reduces the dimension of data while minimizing the loss of information [63]. For the expert opinion dataset (EDS), the authors of the current study selected the soil and terrain indicators for the EDS considering the arid climatic characteristics coupled with low annual perception rate and the associated pedogenic soil processes which plays a major role in forming of soil properties, the gathered knowledge of the evaluated area about the cost of sampling, environmental functions, management practices, vulnerability to productivity, and previous literature recommendations [5, 25, 64–66]. The EDS includes soil electrical conductivity (EC), soil organic matter content (OM), calcium carbonate content (CaCO3), available water holding capacity (AWHC), soil packing density (PD), and topographic witness index (TWI).The soil quality indicators were rescaled to a Z-score with standard normal distribution where they are centered around a mean of zero (0) and standard deviation of 1 as the indicators are of different measurement units and scales. TheZ-score standardization ensures that the indicators have an equal order of magnitude [67].


### 2.4. Minimum Dataset (MDS) Formulation

A minimum dataset (MDS) of soil and terrain indicators was selected from the indicator dataset through multivariate statistical analysis. The minimum dataset selection process includes two steps: representative indicators selection and indicator redundancy reduction. Representative indicators were selected by applying multivariate factor analysis, and the redundant indicators were reduced by vector norm (Norm) analysis.Factor analyses were performed with the covariance matrix of the raw soil and terrain indicators values using principal axis as a factoring method, diagonals equal to 1 as prior communality, and unrotated factor loadings for model parameters. The total dataset includes sand, silt, clay, pH, EC, OM, ESP, CEC, CaCO3, FC, AWHC, BD, PD, Slp, Rlf, and TWI.The raw soil and terrain indicator covariance revealed that the first nine principal components (PCs) account for 0.99% of the total variance of soil and terrain indicators data. The first principal component explains only about 58% of the total variance of soil and terrain indicators, while the second principal components explain about 19% of the total variance. For the other principal components, the contributing indicators were of lowering contribution with a minimum ofr = 0.30.Due to the difference in soil and terrain indicator units of measurements and scales and other than the first and second principal components, the remaining principal components explain only about 23% of the total variability of soil and terrain indicators. This led to the inadequacy of the covariance analyses, and therefore, factor analyses were recomputed based on the correlation matrix of the standardized soil and terrain indicator values where each indicator has a mean centered around zero (0) and standard deviation of unit variance with total variance equal to the number of indicators. Varimax rotation was performed to maximize the relationship between factors and indicators as each indicator has either a small or large loading on each factor [68]. Only factors with eigenvalues more than one were retained as they explain variation in the indicators and not the other factors [69]. For each factor based on the previous criteria, the indicators which contribute with highly weighted factor loadings (within 10% of highest factor loading or absolute value > 0.85) were identified and retained for the minimum dataset (MDS).An indicator value for each representative indicator was calculated, if there were more than one indicator for a single factor to reduce redundancy and exclude pseudo-indicator grouping. At first, for each factor, the soil and terrain indicators that have a factor loading more than or equal to 0.5 were grouped into one group. Any indicator that has all factor loads in the different factor that is more or equal to 0.5 was grouped with the smallest correlations with other indicators. If there is a correlation coefficient less than 0.3 between any indicator and the other indicators, the indicator is placed in another separate group.The factor loadings for selecting soil and terrain indicators for the minimum dataset are not enough as they may ignore some important indicators as the eigenvectors do not express the magnitude (norm) of the resulting factor vectors or for the original indicators [70]. Therefore, for each representative soil and terrain indicator, a vector norm value which represents the magnitude of the vector representing the indicator was calculated [21, 70, 71] according to the following equation:(1)Ni=∑j=1kμij2λi,where Ni is the combined factor load of indicator i in all factors with eigenvalues ≥ 1, uij is the factor load of indicator i in factor j, and λj is the eigenvalue of factor j.Secondly, a regression coefficient of determination (R2) was calculated for each pair of the exterior environmental variables with each of soil and terrain indicators as soil quality influenced by the interior properties and the exterior environment variables [72]. The indicator value was calculated by employing a standard transformation function sum to the vector norm of each indicator with the coefficient of determination R2 for that indicator with the exterior environment variables.


### 2.5. Indicator Weighting

For soil quality index (SQI) formulation, soil and terrain indicators of the three datasets (TDS, MDS, and EDS) were weighted. Each soil and terrain indicator of the total dataset (TDS) and expert opinion dataset (EDS) weighted by calculating the ratio of its respective indicator’s communality resulted from factor analysis and the total indicator communality summation for the dataset [36, 73]. For the minimum dataset (MDS), each indicator was weighted by determining the variation of each respective factor, normalized to unity to the common factor variance summation for the minimum dataset [15] as follows:(2)Iwmds=vijmds∑i=1kvjmds,where Iwmds is the weight for the minimum dataset indicator, vijmds is the variance of the indicator I for factor j, and vjmds is the common factor variance of the minimum dataset indicator.


### 2.6. Indicator Scoring Functions

The soil and terrain indicators are of different measurement units and scales, though there is a need for a normalized transformation of soil and terrain indicator measured values to a unitless score ranging from a zero to one scale [15]. The scoring of indicators provides the capability of combining and averaging the scores into a single value nonlimiting to the pertinent soil functions and processes [74] and to capture information that might otherwise go undetected when examining only the observed values [75]. Linear and nonlinear scoring functions were widely used for quality indicator scoring [15, 76–78].The indicator scoring function relies on critical threshold values which are soil property values determining the upper limit (threshold) where the indicator score is at the most preferable level (score = 1), the lower limit (lower threshold) where the indicator score is at the inadmissible level, and baseline values (minimum target thresholds) where the scoring function equals 0.5 at the lower and upper target thresholds for soil property values [79–82]. The critical threshold values for normalizing indicators values (Table 1) were determined based on the measured values of soil and quality indicators, expert knowledge for the natural ecosystem of the evaluated area, and guidance of the previous literature [80–90].Table 1
Soil quality scoring functions (SSFs) and thresholds values for soil and terrain quality indicators.
IndicatorSSFThresholdsBOL1U1SSourcesLUSandSSF5060—3530500.440[89]SiltSSF302510———0.249[89]ClaySSF303015———0.268[89]pHSSF54.59.5—7.006.57.71.301[3, 75, 83, 84, 90, 96, 97]ECa,bSSF50.02—10.251.752.234[84, 90, 96, 97]OMa,bSSF30.221———1.049[36, 86, 96–98]ESPSSF92155———2.50[15, 39, 66, 99–101]CECaSSF302110.5———0.245[36, 82, 89, 96, 102, 103]CaCO3bSSF9183———2.50NE, EOFCSSF3152520———−2.50NE, EOAWHCa,bSSF32104———0.198[80, 85]; NEBDSSF91.002.01.50———−2.62–0.323[104, 105]PDbSSF91.001.501.2———2.50[81, 97, 106, 107]SlpaSSF92.004.003.00———2.50NE, EORlfSSF93105———2.50NE, EOTWIbSSF951510———2.50NE, EOaSelected soil and terrain indicators based on the minimum dataset (MDS) from factor analysis. bSelected soil and terrain indicators chosen by expert opinion based on the natural ecosystem. NE: natural ecosystem, EO: expert opinion. For SSF3: L, lower threshold, at which or below the score is 0; B, baseline, at which score is 0.5; U, upper threshold, at which or above score is 1.0. For SSF5: L, lower threshold, at which or below the score is 0; L1, lower baseline is 0.5 with bell-shaped relationship; O, optimum level, at which score is 1.0; U1, upper baseline is 0.5 with bell-shaped relationship; and U, upper threshold, at which or above score is 0. For SSF9: L, lower threshold, at which or below the score is 1; B, baseline, at which score is 0.5; and U, upper threshold, at which or above score is 0. S, slope at baseline for nonlinear scoring functions.The soil and terrain indicator datasets (TDS, MDS, and EDS) were scored by transforming data values using linear and nonlinear scoring techniques as described below.


#### 2.6.1. Linear Scoring Functions

Three linear scoring functions (LSFs) were identified and evaluated for scoring soil and terrain indicators: Leibig Linear Scoring Function (LLSF), Homothetic Linear Scoring Function (HLSF), and Glover Linear Scoring Function (GLSF).The soil and terrain indicators were scored according to Liebig scoring function [77] into two orders in terms of soil function: ascending order which is more is better and descending order which is “less is better.” In the ascending order “more is better order,” the indicator value was divided by the maximum indicator values where it is equal to one score (equation (3)). In the descending order “less is better,” the minimum indicator value was divided by the indicator value where the indicator minimum value is equal to one score (equation (4)).(3)LLSFYm=xxmax,(4)LLSFYl=xminx,where X is the indicator value and Xmin and Xmax are the minimum and maximum value of each indicator, respectively.The soil and terrain indicators were scored using a homothetic scoring function [15, 24, 91]. Three types of homothetic scoring functions were used for standardizing the quality indicators to a score value ranging between zero and one: S scoring function (HLSF3), reverse scoring function (HLFS9), and parabola scoring function (HLSF5) [73, 76, 79–82]. For the ascending order “more is better,” HLSF3 scoring function was applied (equation (5)), for “less is better,” HLSF9 scoring function was applied (equation (7)), and for “optimum,” where the indicator value scored as “more is better” until a critical threshold value where it is assigned “less is better,” HLSF5 scoring function was applied (equation (6)).The standard scoring functions (SSFs) are defined as follows:(5)HLSF3:fx=0.1,x≤L,0.1+0.9x−LU−L,L≤x≤U,1.0,x≥U,(6)HLSF5:fx=0.1,x<L,x>U,0.1+0.9x−LU−L,L≤x<L1,1−0.9x−LU−L,U1<x≤U,1.0,L1<x<U1,(7)HLSF9:fx=1.0,x≤L,1−0.9x−LU−L,L≤x≤U,0.1,x≥U,where f(x) is the scoring function and x is the soil property value. For SSF3: L- lower threshold and U- upper threshold. For SSF5: L- lower threshold and L1- lower baseline is 0.5 with a bell-shaped relationship; U1- upper baseline is 0.5 with a bell-shaped relationship; and U- upper threshold. For SSF9: L- lower threshold and U- upper threshold.The soil and terrain indicators were scored using Glover linear scoring function [82, 87, 88]. The quality indicators were scored into two orders in terms of soil function: ascending order which is more is better and descending order which is “less is better.” In the ascending order “more is better order,” equation (8) was used for scoring, and for the descending order “less is better,” equation (9) was used, while for the values outside the lower and upper thresholds, the indicator scored as zero. A combination of both “more is better” and “less is better” was used for the “optimum is better” soil function.(8)GLSFYm=x−st−s,(9)GLSFYl=1−x−st−s,where Ym is the linear score for “more is better,” Yl is the linear score for “less is better,” x is the indicator value, and s and t are lower and upper threshold values.


#### 2.6.2. Nonlinear Scoring Function (NLSF)

Two nonlinear scoring functions (NLSFs) were identified and evaluated for scoring soil and terrain indicators: Sigmoid Nonlinear Scoring Function (SNLSF) and Glover Nonlinear Scoring Function (GNLSF).The soil and terrain indicators were scored using sigmoid nonlinear scoring function [88, 92, 93]. The used sigmoid function for indicator scoring was the logistic nonlinear curve equation (equation (10)) where it is characterized as an S-shaped curve [94]. The scores of the sigmoid function are between zero and one.(10)SNSLFY=a1+x/xob,where Y is the nonlinear score of the indicator value, a is the maximum score, x is the indicator value, xo is the mean value of the indicator, and b is the slope of the curve.The soil and terrain indicators were scored using Glover nonlinear scoring function [82, 87, 88] (equation (11)). Three types of GNLSF scoring were applied: the “more is better” scoring curve for positive slopes, “less is better” scoring curve for negative slopes, and a combination of “Optimum” curve defined by the combination of both “more is better” and “less is better” was used for “optimum is better” soil function [81, 82, 95].(11)GNSLFY=11+e−bx−A,where Y is the nonlinear score of the indicator value, x is the indicator value, A is the baseline value where the indicator score equals 0.5, and b is the slope.


### 2.7. Soil Quality Index Development

Three soil quality indexing methods were computed: the Additive Soil Quality Index (SQI-A), Weighted Additive Soil Quality Index (SQI-W), and Nemoro Soil Quality Index (SQI-N). Each soil quality method was applied for each soil quality indicator dataset.In the Additive Soil Quality Index (SQI-A), the scores of the quality indicators were summed and divided by the dataset number of indicators [15, 108].(12)SQI−A=∑i=1nSin,where SQI-A is the soil quality index, Si is the linear or nonlinear indicator score value, and n is the number of indicators of the dataset.For the Weighted Additive Soil Quality Index (SQI-W), the score of each indicator was multiplied by the indicator weight and summed ([64, 78, 83]).(13)SQI−W=∑i=1nWiSi,where SQI-W is the weighted additive soil quality index, Wi is the indicator weight, Si is the linear or nonlinear indicator score value, and n is the number of indicators of the dataset.The Nemoro Soil Quality Index (SQI-N) was calculated according to the indicator’s minimum and average scores [108–110].(14)SQI−N=Siave2+Simin22×n−1n,where SQI-N is the soil quality index, Siave is the indicator’s score average, Simin is the linear or nonlinear indicator’s score minimum, and n is the number of indicators of the dataset.


### 2.8. Soil Quality Grade Classification

Each soil quality index range was classified into grades by applying Jenks’s natural breaks optimization method [111]. Jenks’s method is an iterative clustering method for determining the best clustering of values into different classes by reducing the variance within classes and maximizing the variance between classes [112]. Five grades were determined for each soil quality index: very high (grade I), high (grade II), moderate (grade III), low (grade IV), and very low (grade V). The grades are from the most suitable for plant growth to the most severe for plant growth.


### 2.9. Soil Quality Index Comparison and Sensitivity

Direct comparison analyses of soil quality grades were used to evaluate the indices’ performance by comparing the number of sampling sites where the soil quality index and indicator combinations had the same soil quality grade [36, 96, 108]. Soil quality indexing methods were also evaluated using sensitivity analysis [88, 113] by computing the ratio between the maximum and minimum soil quality index value for each scoring function using each soil and terrain indicators dataset selection methods. The soil quality indexing method with the higher index value of sensitivity is more preferable as this is sensitive to perplexity in management practices. Also, the linear relationship between the soil quality indices was examined by the correlation within indices, and the regression relation within the indicator methods was performed.


### 2.10. Soil Quality Index Spatial Distribution

The spatial distributions of soil quality index grades were developed by applying the Ordinary Kriging (OK) interpolation technique. For each soil quality index, the semivariogram models were generated using the following model parameters: model type, sill, range, and maximum and minimum neighbors. Eleven model types (circular, spherical, tetraspherical, pentaspherical, exponential, Gaussian, rational quadratic, hole effect, K-Bessel, J-Bessel, and stable) were examined for predicting the spatial distribution of soil quality indices. The accuracy of the predictive model is of greater importance as it is shown how well the spatial pattern of variation can be represented [114]. The best fitting semivariogram model with optimized parameters was selected by incorporating the cross-validation technique to obtain high interpolation accuracy. The cross-validation technique ensures the same probability of being validated against each data sample [115]. The output statistical parameters of the cross-validation technique determine the best and highly accurate prediction model for mapping soil quality indices [115].


### 2.11. Statistical Analyses

Soil indicators were examined for outliers using Quantile Range Outliers statistic with (0.01, 0.05, 0.1, and 0.15) probabilities of the lower quantiles. The normalities of all indicators were examined using the Kolmogorov–Smirnov test and visual examination of histograms. The pH, OM, BD, and CEC data were log transformed, EC was inverse transformed, and SAR, CCE, and ESP were square-root transformed to reduce the skew of their distributions. The relationship between soil and terrain indicators was examined by the Pearson correlation coefficient (r). Bartlett’s test of sphericity was conducted to examine whether the correlation matrix of the soil and terrain indicators is an orthogonal matrix or it diverges significantly from the identity matrix and there is a degree of redundancy between the indicators. The Kaiser–Meyer–Olkin Measure of Sampling Adequacy (MSA) was conducted to indicate whether factor analysis is likely to be appropriate and will gain distinct and reliable factors or not appropriate. The results of Bartlett’s test and MSA will verify if factor analyses can interpret the structure of data and reduce the soil and terrain indicators to a minimum dataset. The statistical and multivariate analyses were performed using JMP pro 14 software for Windows [68]. Scoring and indexing were performed using Microsoft Excel software, ver. 2016 and geostatistical interpolation and mapping using ArcGIS Desktop 10.8 software [116].


## 2.1. Field Description of the Study Area

The evaluated area is located in Wadi El Natrun district, El Beheira Governorate, at the western fringes of the Nile Delta between 29°54′00″E–30°20′00″E and latitudes of 30°22′00″N–30°00′00″N (Figure1). The evaluated area covers an area of 1600 km2 (380000 feddan), with a landscape consisting of Wadi terraces and Wadi depressions with longitudinal sand dunes at the southern edge of the area [43]. The annual evaporation rate is about 114.3 mm/y, the average air temperature is 21°C, and the area has a rarely annual precipitation rate with an average of 41.4 mm/y [44] with an erratic distribution. Based on the work in [45], the soil temperature regime is “Thermic” and the moisture regime is “Torric.” Two soil orders are dominant in the area, Entisols and Aridisols. The evaluated area has been reclaimed during the last two decades and is cultivated with diverse kinds of fruits and vegetables [46]. This kind of cultivation is usually associated with an increase of underground water irrigation, and the ongoing climate change put a threat on the soil quality and sustainable agriculture that reduces the economic value of the lands [47, 48]. Additionally, the excessive extraction of groundwater for irrigation resulted in land subsidence, which could also affect the integrity of the soil profiles in the area [49].Figure 1
Location of the study area and soil sampling sites.


## 2.2. Data Collection

### 2.2.1. Indicator Sampling Design

The required samples size needed for soil quality assessment was determined by using the binomial probability algorithm [50]. 397 sampling sites were chosen. A Proportionate Stratified Random (PSR) sampling technique was used as the sampling type for soil and terrain sampling based on land cover and land use distribution in the evaluated area. The sampling locations were georeferenced (Figure 1) using the Global Positioning System (GPS).


### 2.2.2. Sampling and Analyses

The soil sampling was carried out according to the sampling design. However, some sampling sites were displaced by other sampling sites near the designed sampling locations due to accessibility difficulties, and some certain areas were currently restricted for access, also the presence of elongated sand dunes at the south of the study area. These accessibility limitations resulted in ununiform spatial sampling distribution.Disturbed and undisturbed samples were collected to the depth of the root zone (30 cm) and described according to the work in [51]. Soil samples were crumbled gently by hand without root material and were air dried, grounded, and sieved through a 2 mm sieve. Soil chemical analysis, including pH in 1 : 2.5 soil-water suspension, electrical conductivity (EC) in soil paste extract, cation exchange capacity (CEC), exchangeable sodium percentage (ESP), organic matter (OM), and calcium carbonate were determined according to standard methods of [52]. Soil physical analyses, including particle size distribution using the pipette method, bulk density (BD) using a core method, available water holding capacity (AWHC), and field capacity (FC), were performed according to the work in [53]. Soil packing density (PD) was estimated as a composite index of soil bulk density and soil clay content according to the work in [54].Terrain analyses were performed on a shuttle radar topography mission (SRTM) 1-arc-second v.30 image [55] to obtain some surface parametric information (elevation, slope, relief intensity, and topographic witness index) [56]. The DEM image was firstly sink filled using depression filling algorithm [57]. The slope (Slp) and topographic witness index (TWI) were extracted [58, 59]. The extraction of relief intensity (Rlf) was based on the algorithm described in [60].The Normalized Difference Vegetation Index (NDVI) and elevation data were used as exterior environmental variables in the process of representative indicator selection and redundancy reduction. NDVI data were retrieved from the Landsat Data Continuity Mission (LDCM) sensor multiband dataset using the NDVI equation ([61, 62]).All of the datasets were projected into the WGS84-based Universal Transverse Mercator (UTM) orthographic projection coordinate system (EPSG 32636) and resampled to a 30 m spatial resolution.


## 2.2.1. Indicator Sampling Design

The required samples size needed for soil quality assessment was determined by using the binomial probability algorithm [50]. 397 sampling sites were chosen. A Proportionate Stratified Random (PSR) sampling technique was used as the sampling type for soil and terrain sampling based on land cover and land use distribution in the evaluated area. The sampling locations were georeferenced (Figure 1) using the Global Positioning System (GPS).


## 2.2.2. Sampling and Analyses

The soil sampling was carried out according to the sampling design. However, some sampling sites were displaced by other sampling sites near the designed sampling locations due to accessibility difficulties, and some certain areas were currently restricted for access, also the presence of elongated sand dunes at the south of the study area. These accessibility limitations resulted in ununiform spatial sampling distribution.Disturbed and undisturbed samples were collected to the depth of the root zone (30 cm) and described according to the work in [51]. Soil samples were crumbled gently by hand without root material and were air dried, grounded, and sieved through a 2 mm sieve. Soil chemical analysis, including pH in 1 : 2.5 soil-water suspension, electrical conductivity (EC) in soil paste extract, cation exchange capacity (CEC), exchangeable sodium percentage (ESP), organic matter (OM), and calcium carbonate were determined according to standard methods of [52]. Soil physical analyses, including particle size distribution using the pipette method, bulk density (BD) using a core method, available water holding capacity (AWHC), and field capacity (FC), were performed according to the work in [53]. Soil packing density (PD) was estimated as a composite index of soil bulk density and soil clay content according to the work in [54].Terrain analyses were performed on a shuttle radar topography mission (SRTM) 1-arc-second v.30 image [55] to obtain some surface parametric information (elevation, slope, relief intensity, and topographic witness index) [56]. The DEM image was firstly sink filled using depression filling algorithm [57]. The slope (Slp) and topographic witness index (TWI) were extracted [58, 59]. The extraction of relief intensity (Rlf) was based on the algorithm described in [60].The Normalized Difference Vegetation Index (NDVI) and elevation data were used as exterior environmental variables in the process of representative indicator selection and redundancy reduction. NDVI data were retrieved from the Landsat Data Continuity Mission (LDCM) sensor multiband dataset using the NDVI equation ([61, 62]).All of the datasets were projected into the WGS84-based Universal Transverse Mercator (UTM) orthographic projection coordinate system (EPSG 32636) and resampled to a 30 m spatial resolution.


## 2.3. Indicator Selection and Scaling

Three datasets of quality indicators were chosen, the total dataset (TDS), minimum dataset (MDS), and expert opinion dataset (EDS). The total dataset (TDS) includes all soil and terrain indicators for soil quality index development. The soil and terrain indicators for the minimum dataset were determined according to multivariate factor analysis. Factor analysis reduces the dimension of data while minimizing the loss of information [63]. For the expert opinion dataset (EDS), the authors of the current study selected the soil and terrain indicators for the EDS considering the arid climatic characteristics coupled with low annual perception rate and the associated pedogenic soil processes which plays a major role in forming of soil properties, the gathered knowledge of the evaluated area about the cost of sampling, environmental functions, management practices, vulnerability to productivity, and previous literature recommendations [5, 25, 64–66]. The EDS includes soil electrical conductivity (EC), soil organic matter content (OM), calcium carbonate content (CaCO3), available water holding capacity (AWHC), soil packing density (PD), and topographic witness index (TWI).The soil quality indicators were rescaled to a Z-score with standard normal distribution where they are centered around a mean of zero (0) and standard deviation of 1 as the indicators are of different measurement units and scales. TheZ-score standardization ensures that the indicators have an equal order of magnitude [67].


## 2.4. Minimum Dataset (MDS) Formulation

A minimum dataset (MDS) of soil and terrain indicators was selected from the indicator dataset through multivariate statistical analysis. The minimum dataset selection process includes two steps: representative indicators selection and indicator redundancy reduction. Representative indicators were selected by applying multivariate factor analysis, and the redundant indicators were reduced by vector norm (Norm) analysis.Factor analyses were performed with the covariance matrix of the raw soil and terrain indicators values using principal axis as a factoring method, diagonals equal to 1 as prior communality, and unrotated factor loadings for model parameters. The total dataset includes sand, silt, clay, pH, EC, OM, ESP, CEC, CaCO3, FC, AWHC, BD, PD, Slp, Rlf, and TWI.The raw soil and terrain indicator covariance revealed that the first nine principal components (PCs) account for 0.99% of the total variance of soil and terrain indicators data. The first principal component explains only about 58% of the total variance of soil and terrain indicators, while the second principal components explain about 19% of the total variance. For the other principal components, the contributing indicators were of lowering contribution with a minimum ofr = 0.30.Due to the difference in soil and terrain indicator units of measurements and scales and other than the first and second principal components, the remaining principal components explain only about 23% of the total variability of soil and terrain indicators. This led to the inadequacy of the covariance analyses, and therefore, factor analyses were recomputed based on the correlation matrix of the standardized soil and terrain indicator values where each indicator has a mean centered around zero (0) and standard deviation of unit variance with total variance equal to the number of indicators. Varimax rotation was performed to maximize the relationship between factors and indicators as each indicator has either a small or large loading on each factor [68]. Only factors with eigenvalues more than one were retained as they explain variation in the indicators and not the other factors [69]. For each factor based on the previous criteria, the indicators which contribute with highly weighted factor loadings (within 10% of highest factor loading or absolute value > 0.85) were identified and retained for the minimum dataset (MDS).An indicator value for each representative indicator was calculated, if there were more than one indicator for a single factor to reduce redundancy and exclude pseudo-indicator grouping. At first, for each factor, the soil and terrain indicators that have a factor loading more than or equal to 0.5 were grouped into one group. Any indicator that has all factor loads in the different factor that is more or equal to 0.5 was grouped with the smallest correlations with other indicators. If there is a correlation coefficient less than 0.3 between any indicator and the other indicators, the indicator is placed in another separate group.The factor loadings for selecting soil and terrain indicators for the minimum dataset are not enough as they may ignore some important indicators as the eigenvectors do not express the magnitude (norm) of the resulting factor vectors or for the original indicators [70]. Therefore, for each representative soil and terrain indicator, a vector norm value which represents the magnitude of the vector representing the indicator was calculated [21, 70, 71] according to the following equation:(1)Ni=∑j=1kμij2λi,where Ni is the combined factor load of indicator i in all factors with eigenvalues ≥ 1, uij is the factor load of indicator i in factor j, and λj is the eigenvalue of factor j.Secondly, a regression coefficient of determination (R2) was calculated for each pair of the exterior environmental variables with each of soil and terrain indicators as soil quality influenced by the interior properties and the exterior environment variables [72]. The indicator value was calculated by employing a standard transformation function sum to the vector norm of each indicator with the coefficient of determination R2 for that indicator with the exterior environment variables.


## 2.5. Indicator Weighting

For soil quality index (SQI) formulation, soil and terrain indicators of the three datasets (TDS, MDS, and EDS) were weighted. Each soil and terrain indicator of the total dataset (TDS) and expert opinion dataset (EDS) weighted by calculating the ratio of its respective indicator’s communality resulted from factor analysis and the total indicator communality summation for the dataset [36, 73]. For the minimum dataset (MDS), each indicator was weighted by determining the variation of each respective factor, normalized to unity to the common factor variance summation for the minimum dataset [15] as follows:(2)Iwmds=vijmds∑i=1kvjmds,where Iwmds is the weight for the minimum dataset indicator, vijmds is the variance of the indicator I for factor j, and vjmds is the common factor variance of the minimum dataset indicator.


## 2.6. Indicator Scoring Functions

The soil and terrain indicators are of different measurement units and scales, though there is a need for a normalized transformation of soil and terrain indicator measured values to a unitless score ranging from a zero to one scale [15]. The scoring of indicators provides the capability of combining and averaging the scores into a single value nonlimiting to the pertinent soil functions and processes [74] and to capture information that might otherwise go undetected when examining only the observed values [75]. Linear and nonlinear scoring functions were widely used for quality indicator scoring [15, 76–78].The indicator scoring function relies on critical threshold values which are soil property values determining the upper limit (threshold) where the indicator score is at the most preferable level (score = 1), the lower limit (lower threshold) where the indicator score is at the inadmissible level, and baseline values (minimum target thresholds) where the scoring function equals 0.5 at the lower and upper target thresholds for soil property values [79–82]. The critical threshold values for normalizing indicators values (Table 1) were determined based on the measured values of soil and quality indicators, expert knowledge for the natural ecosystem of the evaluated area, and guidance of the previous literature [80–90].Table 1
Soil quality scoring functions (SSFs) and thresholds values for soil and terrain quality indicators.
IndicatorSSFThresholdsBOL1U1SSourcesLUSandSSF5060—3530500.440[89]SiltSSF302510———0.249[89]ClaySSF303015———0.268[89]pHSSF54.59.5—7.006.57.71.301[3, 75, 83, 84, 90, 96, 97]ECa,bSSF50.02—10.251.752.234[84, 90, 96, 97]OMa,bSSF30.221———1.049[36, 86, 96–98]ESPSSF92155———2.50[15, 39, 66, 99–101]CECaSSF302110.5———0.245[36, 82, 89, 96, 102, 103]CaCO3bSSF9183———2.50NE, EOFCSSF3152520———−2.50NE, EOAWHCa,bSSF32104———0.198[80, 85]; NEBDSSF91.002.01.50———−2.62–0.323[104, 105]PDbSSF91.001.501.2———2.50[81, 97, 106, 107]SlpaSSF92.004.003.00———2.50NE, EORlfSSF93105———2.50NE, EOTWIbSSF951510———2.50NE, EOaSelected soil and terrain indicators based on the minimum dataset (MDS) from factor analysis. bSelected soil and terrain indicators chosen by expert opinion based on the natural ecosystem. NE: natural ecosystem, EO: expert opinion. For SSF3: L, lower threshold, at which or below the score is 0; B, baseline, at which score is 0.5; U, upper threshold, at which or above score is 1.0. For SSF5: L, lower threshold, at which or below the score is 0; L1, lower baseline is 0.5 with bell-shaped relationship; O, optimum level, at which score is 1.0; U1, upper baseline is 0.5 with bell-shaped relationship; and U, upper threshold, at which or above score is 0. For SSF9: L, lower threshold, at which or below the score is 1; B, baseline, at which score is 0.5; and U, upper threshold, at which or above score is 0. S, slope at baseline for nonlinear scoring functions.The soil and terrain indicator datasets (TDS, MDS, and EDS) were scored by transforming data values using linear and nonlinear scoring techniques as described below.


### 2.6.1. Linear Scoring Functions

Three linear scoring functions (LSFs) were identified and evaluated for scoring soil and terrain indicators: Leibig Linear Scoring Function (LLSF), Homothetic Linear Scoring Function (HLSF), and Glover Linear Scoring Function (GLSF).The soil and terrain indicators were scored according to Liebig scoring function [77] into two orders in terms of soil function: ascending order which is more is better and descending order which is “less is better.” In the ascending order “more is better order,” the indicator value was divided by the maximum indicator values where it is equal to one score (equation (3)). In the descending order “less is better,” the minimum indicator value was divided by the indicator value where the indicator minimum value is equal to one score (equation (4)).(3)LLSFYm=xxmax,(4)LLSFYl=xminx,where X is the indicator value and Xmin and Xmax are the minimum and maximum value of each indicator, respectively.The soil and terrain indicators were scored using a homothetic scoring function [15, 24, 91]. Three types of homothetic scoring functions were used for standardizing the quality indicators to a score value ranging between zero and one: S scoring function (HLSF3), reverse scoring function (HLFS9), and parabola scoring function (HLSF5) [73, 76, 79–82]. For the ascending order “more is better,” HLSF3 scoring function was applied (equation (5)), for “less is better,” HLSF9 scoring function was applied (equation (7)), and for “optimum,” where the indicator value scored as “more is better” until a critical threshold value where it is assigned “less is better,” HLSF5 scoring function was applied (equation (6)).The standard scoring functions (SSFs) are defined as follows:(5)HLSF3:fx=0.1,x≤L,0.1+0.9x−LU−L,L≤x≤U,1.0,x≥U,(6)HLSF5:fx=0.1,x<L,x>U,0.1+0.9x−LU−L,L≤x<L1,1−0.9x−LU−L,U1<x≤U,1.0,L1<x<U1,(7)HLSF9:fx=1.0,x≤L,1−0.9x−LU−L,L≤x≤U,0.1,x≥U,where f(x) is the scoring function and x is the soil property value. For SSF3: L- lower threshold and U- upper threshold. For SSF5: L- lower threshold and L1- lower baseline is 0.5 with a bell-shaped relationship; U1- upper baseline is 0.5 with a bell-shaped relationship; and U- upper threshold. For SSF9: L- lower threshold and U- upper threshold.The soil and terrain indicators were scored using Glover linear scoring function [82, 87, 88]. The quality indicators were scored into two orders in terms of soil function: ascending order which is more is better and descending order which is “less is better.” In the ascending order “more is better order,” equation (8) was used for scoring, and for the descending order “less is better,” equation (9) was used, while for the values outside the lower and upper thresholds, the indicator scored as zero. A combination of both “more is better” and “less is better” was used for the “optimum is better” soil function.(8)GLSFYm=x−st−s,(9)GLSFYl=1−x−st−s,where Ym is the linear score for “more is better,” Yl is the linear score for “less is better,” x is the indicator value, and s and t are lower and upper threshold values.


### 2.6.2. Nonlinear Scoring Function (NLSF)

Two nonlinear scoring functions (NLSFs) were identified and evaluated for scoring soil and terrain indicators: Sigmoid Nonlinear Scoring Function (SNLSF) and Glover Nonlinear Scoring Function (GNLSF).The soil and terrain indicators were scored using sigmoid nonlinear scoring function [88, 92, 93]. The used sigmoid function for indicator scoring was the logistic nonlinear curve equation (equation (10)) where it is characterized as an S-shaped curve [94]. The scores of the sigmoid function are between zero and one.(10)SNSLFY=a1+x/xob,where Y is the nonlinear score of the indicator value, a is the maximum score, x is the indicator value, xo is the mean value of the indicator, and b is the slope of the curve.The soil and terrain indicators were scored using Glover nonlinear scoring function [82, 87, 88] (equation (11)). Three types of GNLSF scoring were applied: the “more is better” scoring curve for positive slopes, “less is better” scoring curve for negative slopes, and a combination of “Optimum” curve defined by the combination of both “more is better” and “less is better” was used for “optimum is better” soil function [81, 82, 95].(11)GNSLFY=11+e−bx−A,where Y is the nonlinear score of the indicator value, x is the indicator value, A is the baseline value where the indicator score equals 0.5, and b is the slope.


## 2.6.1. Linear Scoring Functions

Three linear scoring functions (LSFs) were identified and evaluated for scoring soil and terrain indicators: Leibig Linear Scoring Function (LLSF), Homothetic Linear Scoring Function (HLSF), and Glover Linear Scoring Function (GLSF).The soil and terrain indicators were scored according to Liebig scoring function [77] into two orders in terms of soil function: ascending order which is more is better and descending order which is “less is better.” In the ascending order “more is better order,” the indicator value was divided by the maximum indicator values where it is equal to one score (equation (3)). In the descending order “less is better,” the minimum indicator value was divided by the indicator value where the indicator minimum value is equal to one score (equation (4)).(3)LLSFYm=xxmax,(4)LLSFYl=xminx,where X is the indicator value and Xmin and Xmax are the minimum and maximum value of each indicator, respectively.The soil and terrain indicators were scored using a homothetic scoring function [15, 24, 91]. Three types of homothetic scoring functions were used for standardizing the quality indicators to a score value ranging between zero and one: S scoring function (HLSF3), reverse scoring function (HLFS9), and parabola scoring function (HLSF5) [73, 76, 79–82]. For the ascending order “more is better,” HLSF3 scoring function was applied (equation (5)), for “less is better,” HLSF9 scoring function was applied (equation (7)), and for “optimum,” where the indicator value scored as “more is better” until a critical threshold value where it is assigned “less is better,” HLSF5 scoring function was applied (equation (6)).The standard scoring functions (SSFs) are defined as follows:(5)HLSF3:fx=0.1,x≤L,0.1+0.9x−LU−L,L≤x≤U,1.0,x≥U,(6)HLSF5:fx=0.1,x<L,x>U,0.1+0.9x−LU−L,L≤x<L1,1−0.9x−LU−L,U1<x≤U,1.0,L1<x<U1,(7)HLSF9:fx=1.0,x≤L,1−0.9x−LU−L,L≤x≤U,0.1,x≥U,where f(x) is the scoring function and x is the soil property value. For SSF3: L- lower threshold and U- upper threshold. For SSF5: L- lower threshold and L1- lower baseline is 0.5 with a bell-shaped relationship; U1- upper baseline is 0.5 with a bell-shaped relationship; and U- upper threshold. For SSF9: L- lower threshold and U- upper threshold.The soil and terrain indicators were scored using Glover linear scoring function [82, 87, 88]. The quality indicators were scored into two orders in terms of soil function: ascending order which is more is better and descending order which is “less is better.” In the ascending order “more is better order,” equation (8) was used for scoring, and for the descending order “less is better,” equation (9) was used, while for the values outside the lower and upper thresholds, the indicator scored as zero. A combination of both “more is better” and “less is better” was used for the “optimum is better” soil function.(8)GLSFYm=x−st−s,(9)GLSFYl=1−x−st−s,where Ym is the linear score for “more is better,” Yl is the linear score for “less is better,” x is the indicator value, and s and t are lower and upper threshold values.


## 2.6.2. Nonlinear Scoring Function (NLSF)

Two nonlinear scoring functions (NLSFs) were identified and evaluated for scoring soil and terrain indicators: Sigmoid Nonlinear Scoring Function (SNLSF) and Glover Nonlinear Scoring Function (GNLSF).The soil and terrain indicators were scored using sigmoid nonlinear scoring function [88, 92, 93]. The used sigmoid function for indicator scoring was the logistic nonlinear curve equation (equation (10)) where it is characterized as an S-shaped curve [94]. The scores of the sigmoid function are between zero and one.(10)SNSLFY=a1+x/xob,where Y is the nonlinear score of the indicator value, a is the maximum score, x is the indicator value, xo is the mean value of the indicator, and b is the slope of the curve.The soil and terrain indicators were scored using Glover nonlinear scoring function [82, 87, 88] (equation (11)). Three types of GNLSF scoring were applied: the “more is better” scoring curve for positive slopes, “less is better” scoring curve for negative slopes, and a combination of “Optimum” curve defined by the combination of both “more is better” and “less is better” was used for “optimum is better” soil function [81, 82, 95].(11)GNSLFY=11+e−bx−A,where Y is the nonlinear score of the indicator value, x is the indicator value, A is the baseline value where the indicator score equals 0.5, and b is the slope.


## 2.7. Soil Quality Index Development

Three soil quality indexing methods were computed: the Additive Soil Quality Index (SQI-A), Weighted Additive Soil Quality Index (SQI-W), and Nemoro Soil Quality Index (SQI-N). Each soil quality method was applied for each soil quality indicator dataset.In the Additive Soil Quality Index (SQI-A), the scores of the quality indicators were summed and divided by the dataset number of indicators [15, 108].(12)SQI−A=∑i=1nSin,where SQI-A is the soil quality index, Si is the linear or nonlinear indicator score value, and n is the number of indicators of the dataset.For the Weighted Additive Soil Quality Index (SQI-W), the score of each indicator was multiplied by the indicator weight and summed ([64, 78, 83]).(13)SQI−W=∑i=1nWiSi,where SQI-W is the weighted additive soil quality index, Wi is the indicator weight, Si is the linear or nonlinear indicator score value, and n is the number of indicators of the dataset.The Nemoro Soil Quality Index (SQI-N) was calculated according to the indicator’s minimum and average scores [108–110].(14)SQI−N=Siave2+Simin22×n−1n,where SQI-N is the soil quality index, Siave is the indicator’s score average, Simin is the linear or nonlinear indicator’s score minimum, and n is the number of indicators of the dataset.


## 2.8. Soil Quality Grade Classification

Each soil quality index range was classified into grades by applying Jenks’s natural breaks optimization method [111]. Jenks’s method is an iterative clustering method for determining the best clustering of values into different classes by reducing the variance within classes and maximizing the variance between classes [112]. Five grades were determined for each soil quality index: very high (grade I), high (grade II), moderate (grade III), low (grade IV), and very low (grade V). The grades are from the most suitable for plant growth to the most severe for plant growth.


## 2.9. Soil Quality Index Comparison and Sensitivity

Direct comparison analyses of soil quality grades were used to evaluate the indices’ performance by comparing the number of sampling sites where the soil quality index and indicator combinations had the same soil quality grade [36, 96, 108]. Soil quality indexing methods were also evaluated using sensitivity analysis [88, 113] by computing the ratio between the maximum and minimum soil quality index value for each scoring function using each soil and terrain indicators dataset selection methods. The soil quality indexing method with the higher index value of sensitivity is more preferable as this is sensitive to perplexity in management practices. Also, the linear relationship between the soil quality indices was examined by the correlation within indices, and the regression relation within the indicator methods was performed.


## 2.10. Soil Quality Index Spatial Distribution

The spatial distributions of soil quality index grades were developed by applying the Ordinary Kriging (OK) interpolation technique. For each soil quality index, the semivariogram models were generated using the following model parameters: model type, sill, range, and maximum and minimum neighbors. Eleven model types (circular, spherical, tetraspherical, pentaspherical, exponential, Gaussian, rational quadratic, hole effect, K-Bessel, J-Bessel, and stable) were examined for predicting the spatial distribution of soil quality indices. The accuracy of the predictive model is of greater importance as it is shown how well the spatial pattern of variation can be represented [114]. The best fitting semivariogram model with optimized parameters was selected by incorporating the cross-validation technique to obtain high interpolation accuracy. The cross-validation technique ensures the same probability of being validated against each data sample [115]. The output statistical parameters of the cross-validation technique determine the best and highly accurate prediction model for mapping soil quality indices [115].


## 2.11. Statistical Analyses

Soil indicators were examined for outliers using Quantile Range Outliers statistic with (0.01, 0.05, 0.1, and 0.15) probabilities of the lower quantiles. The normalities of all indicators were examined using the Kolmogorov–Smirnov test and visual examination of histograms. The pH, OM, BD, and CEC data were log transformed, EC was inverse transformed, and SAR, CCE, and ESP were square-root transformed to reduce the skew of their distributions. The relationship between soil and terrain indicators was examined by the Pearson correlation coefficient (r). Bartlett’s test of sphericity was conducted to examine whether the correlation matrix of the soil and terrain indicators is an orthogonal matrix or it diverges significantly from the identity matrix and there is a degree of redundancy between the indicators. The Kaiser–Meyer–Olkin Measure of Sampling Adequacy (MSA) was conducted to indicate whether factor analysis is likely to be appropriate and will gain distinct and reliable factors or not appropriate. The results of Bartlett’s test and MSA will verify if factor analyses can interpret the structure of data and reduce the soil and terrain indicators to a minimum dataset. The statistical and multivariate analyses were performed using JMP pro 14 software for Windows [68]. Scoring and indexing were performed using Microsoft Excel software, ver. 2016 and geostatistical interpolation and mapping using ArcGIS Desktop 10.8 software [116].


## 3. Results and Discussion

### 3.1. Soil and Terrain Indicator Characteristics

The soil and terrain indicators related to soil quality were subjected to descriptive analyses (Table2). The quantile range outlier’s statistic indicated no outliers existed in soil and terrain indicator data. The elevation of the sampling sites ranged between −19 and 197 m above mean sea level. The Rlf, sand, Slp, CaCO3, and clay indicators showed high data dispersion indicating that the spatial variations of these soil indicators are large, while the other soil and terrain indicators are of weak spatial variation as they have low data dispersion.Table 2
Descriptive statistics of soil and terrain indicators.
IndicatorMeanSEMinQ1MedianQ3MaxSand81.090.3965.0076.0080.0085.0099.87Silt6.160.170.095.007.5410.0017.79Clay8.510.280.048.0012.0015.0026.00pH7.380.027.007.167.247.508.20EC2.630.100.052.022.963.8613.14OM0.220.010.010.140.210.371.50ESP6.430.161.004.516.929.9813.46CEC2.130.100.021.462.323.1611.71CaCO31.950.130.041.442.634.2017.22FC10.950.232.899.5412.3915.0224.63AWHC5.400.062.864.895.606.358.61BD1.590.001.411.591.611.621.63PD1.690.001.421.671.721.751.84Slp5.900.220.544.096.329.8924.11Rlf14.940.283.0012.0015.0019.0044.00TWI7.790.095.526.827.248.3311.70Kaiser–Meyer–Olkin measure of sampling adequacyMSA0.74Bartlett’s test of sphericityApprox. chi-square16428.05Df120Sig.0.00Mean is represented by geometric mean, SE is the standard error, Min is the minimum value, Max is the maximum value,Q1 and Q3 are the lower and upper quartiles, respectively.The correlation analysis of soil physical and chemical indicators and terrain indicators revealed that there is a significant correlation in 35 of 136 soil and terrain indicator correlation pairs at a significant level ofp<0.05 (Table 3). The highest positive correlations were between soil clay content against soil field capacity, available water holding capacity, and soil packing density where the correlation coefficient is greater than 0.90. Also, the soil field capacity indicated a high positive correlation coefficient (r > 0.95) with available water holding capacity and soil packing density (r > 0.94). The highest negative correlation coefficients were obtained between sand and clay (r > 0.93), soil field capacity (r > 0.96), available water holding capacity (r > 0.99), and soil packing density (r > 0.93). It is also noted that soil organic matter indicator (OM) has a moderate positive relationship with CEC (0.45) and a moderate negative relationship with BD (−0.54) while it has no significant relationship with the soil salinity indicator (EC).Table 3
Correlation matrix for measured soil and terrain attributes (n = 397).
SandSiltClaypHECOMESPCECCaCO3FCAWHCBDPDSlpRlfSand1.00Silt−0.791.00Clay−0.930.511.00pH0.70−0.44−0.721.00EC−0.180.130.18−0.271.00OM0.08−0.09−0.060.31−0.031.00ESP−0.050.000.06−0.210.35−0.201.00CEC−0.040.030.040.270.020.45−0.081.00CaCO30.07−0.05−0.060.11−0.040.15−0.070.491.00FC−0.960.590.99−0.690.180.010.030.07−0.061.00AWHC−0.990.800.91−0.640.180.020.000.09−0.060.951.00BD−0.830.590.81−0.780.16−0.540.17−0.16−0.100.780.741.00PD−0.930.570.96−0.790.18−0.310.12−0.06−0.090.940.870.941.00Slp0.03−0.04−0.020.040.00−0.01−0.06−0.03−0.01−0.02−0.03−0.03−0.031.00Rlf0.04−0.01−0.040.10−0.020.08−0.110.080.06−0.03−0.02−0.09−0.060.511.00TWI0.06−0.13−0.010.000.06−0.130.09−0.14−0.12−0.03−0.080.040.02−0.14−0.04EC, electrical conductivity (dS m1); SOM, soil organic matter (%); ESP, exchangeable sodium percent (%); CEC, cation exchange capacity (cmol kg1); CaCO3, calcium carbonates content (%); FC, field capacity water content (%); AWHC, available water holding capacity (%); BD, soil bulk density; PD, soil packing density; Slp, slope (%); Rlf, relief (m), TWI, topographic wetness index.Bartlett’s test of sphericity revealed that the soil and terrain indicator correlation coefficients are significantly different than the identity correlation matrix (<0.05) and the indicators are not perfectly uncorrelated. The Kaiser–Meyer–Olkin Measure of Sampling Adequacy (MSA) revealed that the MSA value is above 0.6 and acceptable. This result of MSA indicates that there is a proportion of variance in soil and terrain indicators that might be caused by underlying factors. The Kaiser–Meyer–Olkin and Bartlett’s test show that there is a high scope for data reduction with factor analyzing where the pattern of data structure in indicators can be interpreted.


### 3.2. Minimum Dataset Indicators

The analyses of the relationship between soil and terrain indicators revealed that there are noted redundancy and collinearity among indicators and there is a need to establish a minimum dataset for the process of soil quality assessment by eliminating some indicators performing factor analyses.The results of factor analyses on standardized soil and terrain indicators values revealed that the first five factors accounted for 78.76% of indicator variability with eigenvalue greater than 1 (Table4). Thus, the first five factor analyses were retained and the other factors with eigenvalue less than 1 were excluded as the indicators may explain more variance than those excluded factors [117].Table 4
Eigenvalues, proportions and cumulative variance explained by factor analysis using a correlation matrix of standardized soil and terrain indicator data.
FactorEigenvalueDifferenceProportionCumulative16.726.5540.9440.9422.091.6510.3451.2831.521.539.5860.8641.281.459.0869.9451.001.418.8278.76The retained factors explain more than 97% of the variance for sand and AWHC in indicators (Table5); more than 91% for PD, FC, BD, and clay indicators; more than 71% for OM, CaCO3, Slp, CEC, pH, Rlf, and EC indicators; and more than 55% for ESP and silt indicators. The TWI indicator was the least important indicator due to its lowest communality estimate (less than 28% variance) as the higher communality estimates procure more preference upon the lower ones.Table 5
Proportion of variance using varimax rotation and communality estimates for soil and terrain indicators for each of the retained factors.
IndicatorFactor 1Factor 2Factor 3Factor 4Factor 5Communality estimatesSand−0.99−0.020.01−0.01−0.030.99Silt0.740.020.000.05−0.050.55Clay0.950.01−0.01−0.010.070.91pH−0.740.350.040.10−0.220.74EC0.170.130.05−0.070.810.71OM−0.090.920.030.09−0.080.87ESP0.00−0.23−0.09−0.020.790.68CEC0.030.55−0.030.670.070.75CaCO3-0.080.03−0.050.880.040.78FC0.970.09−0.01-0.010.050.96AWHC0.980.130.000.020.000.97BD0.85−0.47-0.040.000.100.95PD0.95−0.23−0.02−0.010.090.96Slp−0.02−0.070.880.02−0.020.78Rlf−0.030.090.850.06−0.020.73TWI−0.070.01−0.17−0.450.210.28The representative soil and terrain indicators for soil quality were chosen based on their relationship with the retained factors as the eigenvalues were the criteria magnitude for the interpretation (Table5). In the first factor, eight quality indicator factor loads are more than 0.7 which are sand, silt, clay, pH, FC, AWHC, BD, and PD. In the second factor, only the OM indicator is of high weight variable accounting for 0.92 of variation. In the third factor, slp and rlf indicators’ factor loads are of more than 0.85. CEC and CaCO3 indicators are more than 0.67 in the fourth factor. The fifth factor contains only EC and ESP indicator factor loads of more than 0.79. The soil and terrain quality indicators were grouped according to that criterion into five groups (Table 6).Table 6
Included soil and terrain indicators in the minimum dataset from representative indicators.
GroupIndicatorVector normCoefficient of determinationStandard normal transformationIndicator valueIncludedNDVIElevNormNDVIElev1Sand2.570.180.141.000.801.002.80No1Silt1.910.100.100.740.430.741.91No1Clay2.470.170.110.960.740.962.66No1pH2.010.020.010.780.110.781.67No1FC2.520.200.140.980.900.982.86No1AWHC2.540.230.170.991.000.992.97Yes1BD2.300.020.020.890.090.891.87No1PD2.480.090.070.960.390.962.32No2OM1.360.090.020.530.380.531.43Yes3Slp1.090.000.010.420.020.420.86Yes3Rlf1.060.000.010.410.010.410.83No4CEC1.100.000.010.430.020.430.87Yes4CaCO31.020.000.000.400.010.400.80No5EC0.950.010.000.370.030.370.77Yes5ESP0.860.010.060.330.020.330.69NoThe first factor explained about 41% of the total variance with high positive loadings from AWHC, FC, clay, PD, and BD indicators (>0.85) and high negative loading from the sand indicator (−0.99). The moderate positive loadings were from silt (0.74) and negative loading from pH (−0.74) indicators. The significant correlation coefficients among sand, silt, clay, pH, FC, AWHC, BD, and PD indicators influenced the factors’ loading values. The first factor expresses the influence of soil texture on hydraulic soil properties as it clarifies the soil moisture and water flow properties.The second factor explained about 10% of the total variance with high positive loading from the OM indicator (0.92) and moderate positive loading from CEC (0.55). The second factor had also moderate negative loading from the BD indicator (−0.47) and moderate positive loading from the pH indicator (0.35). The second factor expresses the effect of soil organic matter content on bulk density and soil cation exchange capacity.The third factor explained about 9.6% of the total variance with high positive loading from Slp (0.88) and Rlf (0.85) indicators. This factor is a proxy of surface as it clearly shows the effect of terrain on soil properties as they significantly correlated.The fourth factor explained about 8% of the total variance with high positive loading from the CaCO3 indicator (0.88) and moderate positive loading from CEC (0.67) indicators. Also, the fourth factor has negative moderate loading from the TWI indicator (−0.45) resulting from significant correlations among CaCO3, CEC, and TWI. The fourth factor expresses the interaction between CaCO3, CEC, and soil wetness as where wetness occurs, the calcium salts readily dissolved and interfere with soil exchangeable cations [118, 119] and the presence of calcium carbonates is associated with cascading changes in soil biogeochemistry [120].The fifth factor explained about 8.8% of the total variance with high positive loading from the EC indicator (0.81) and ESP (0.79). This factor represents salinity characteristics.The Norm analyses (Table6) revealed that, from eight quality indicators in the first group, only the AWHC indicator is with the highest Norm value. The OM indicator is the only indicator in the second group, while Slp, CEC, and EC indicators are the highest in the third, fourth, and fifth group, respectively. According to these results, the minimum dataset indicators include five soil and terrain quality indicators which are AWHC, OM, Slp, CEC, and EC indicators.


### 3.3. Soil Quality Assessment

The soil quality indices were calculated and ranked for each method. The classification criteria for the different calculated soil quality indices based on the different scoring techniques show variation between quality methods (A, W, and N), but there are slight similarities between each quality scoring technique for the different quality calculation methods. The threshold grade values of the scoring methods for each soil quality method were nearly the same threshold range values although the threshold grades values for the different scoring methods were not consistent.The cross validation for comparing the different selected interpolation ordinary kriging models identified the exponential predictive model as the best accurate model for representing the spatial distribution of soil quality indices other than the other predictive models (data not shown). The distribution of soil quality indices (Table7) shows a degree of pattern similarity between soil quality indices across the different quality grouping types with indicators scoring functions.Table 7
Distribution of index area grades across different soil quality index methods.
Indicator datasetsSQI modelStatisticSoil quality gradeI (very high)II (high)III (moderate)IV (low)V (very low)TDSSQI-AMin0.0013.8129.7312.340.12Max3.5231.0044.6935.823.53Mean1.7723.7039.6923.611.52SQI-WMin2.8018.0037.773.310.00Max5.7542.1949.7218.870.44Mean4.3135.1342.008.770.09SQI-NMin0.0016.8936.2612.110.11Max2.8425.3251.3728.081.52Mean0.8422.0246.6419.880.91MDSSQI-AMin0.0013.9319.480.170.00Max6.0863.0046.8933.350.17Mean2.4445.7433.828.260.03SQI-WMin0.0112.0429.750.400.00Max3.1258.4646.6236.550.08Mean1.0041.2238.099.960.02SQI-NMin0.001.9418.280.000.00Max3.9458.3247.2767.422.65Mean1.5337.4034.9915.830.55EDSSQI-AMin0.0013.2031.6813.110.00Max5.8828.8446.5237.504.04Mean1.5721.7840.5525.101.30SQI-WMin0.0012.5437.3812.640.00Max6.7828.3347.8128.816.08Mean1.7820.3742.9722.612.56SQI-NMin0.009.1438.7816.330.00Max5.6731.1553.0936.113.12Mean1.4322.6343.1222.410.70Soil quality grade III (moderate) is mainly the largest area percentage proportion in most soil quality indexing methods and scoring functions across the three soil quality datasets groups (TDS, MDS, and EDS) followed by grade II (High), grade IV (low), grade I (very high), and grade V (very low), respectively. A little part of the study area is for grades I (very high) and V (very low) with area percentage coverage. Soil quality grade III (moderate) on average accounts for 40% of the study area, and soil quality grade II (high) accounts for 30% of the study area, while the rest of soil quality grades (I, IV, and V) accounts for 20% of the study area. The disparities between the spatial distribution of soil qualities across indexing methods and scoring functions are of no large distribution disparity.


### 3.4. Comparison of Indices

Sensitivity analyses for soil quality indices showed that the soil quality grouping type of the indicator dataset, quality indexing method, and indicators scoring function are impacting the soil quality values in the study area (Table8). The LLSF linear scoring function is larger sensitive to the variation in soil quality compared with other sensitivities for all indicator grouping types, soil quality methods, and scoring functions specially with the minimum dataset grouping type while SNLSF nonlinear scoring function is the less sensitive to soil quality variation.Table 8
Sensitivity analyses between indicator datasets, models, and scoring functions.
Indicator datasetsSQI modelLSFNLSFLLSFHLSFGLSFSNLSFGNLSFTDSSQI-A2.481.771.621.392.14SQI-W3.683.191.851.532.68SQI-N2.521.891.821.692.14MDSSQI-A4.663.133.312.522.63SQI-W4.874.433.492.543.11SQI-N4.443.433.233.183.52EDSSQI-A3.612.582.212.023.38SQI-W3.943.722.462.422.83SQI-N3.383.063.282.193.09Sensitivity analysis showed that both indexing and dataset selection methods influenced the SQI values for the various soil quality values. Thus, the analyses emphasize that the indicator dataset grouping type, the indexing method, and scoring function could assist in the evaluation of the soil condition status in the study area. The Linear Scoring Functions (LSFs) disclosed more soil quality variations further than Nonlinear Scoring Functions (NLSFs). The resulted sensitivity analyses agree with the results in [89] where the Linear Scoring functions (LSFs) are more favorable than Nonlinear Scoring Functions (NLSFs) in evaluating soil quality indices. It is in contrast with the work in [15, 88] where it was reported that Nonlinear Scoring Functions (NLSFs) are favorable in studies with one linear and one nonlinear scoring functions, and their sensitivity was not evaluated [89].The linear relationships between the indicator datasets with different scoring functions and indices models (Table9) showed a significant correlation (P<0.05). The values of the linear relationships between MDS and TDS datasets are higher than the values of the linear relationships between EDS and TDS datasets across linear and nonlinear scoring functions and index models. The values of the linear relationships between TDS and MDS datasets for the weighted additive quality index are higher than the values of the linear relationship for additive and Nemoro indexing methods for all the linear and nonlinear scoring functions. In contrast, the linear relationships between the TDS and EDS datasets for the weighted additive quality index are mostly lower than the linear relationship values for additive and Nemoro indexing methods except for LLSF and HLSF linear scoring functions. The additive and weighted additive quality models’ values are quite similar for the scoring function type between the datasets.Table 9
The linear relationships between indicator datasets with different scoring functions.
SQIScoring methodSQI modelRelationshipR2SQIScoring methodSQI modelRelationshipR2TDS-MDSLLSFSQI-A0.4586x + 0.31320.54TDS-EDSLLSFSQI-A0.4634x + 0.23780.30SQI-W0.895x + 0.01080.90SQI-W1.0574x − 0.18690.64SQI-N0.5027x + 0.21110.51SQI-N0.4161x + 0.18650.23HLSFSQI-A0.3592x + 0.28880.46HLSFSQI-A0.4433x + 0.20450.43SQI-W0.61x + 0.06610.62SQI-W0.6876x − 0.01290.50SQI-N0.2621x + 0.19090.49SQI-N0.5026x + 0.13930.42GLSFSQI-A0.3725x + 0.24090.42GLSFSQI-A0.4703x + 0.150.39SQI-W0.4989x + 0.06650.51SQI-W0.5191x + 0.01430.31SQI-N0.4239x + 0.1610.41SQI-N0.5284x + 0.09970.39SNLSFSQI-A0.3487x + 0.3210.42SNLSFSQI-A0.3989x + 0.29960.50SQI-W0.4977x + 0.14890.48SQI-W0.4134x + 0.19260.32SQI-N0.3606x + 0.23190.39SQI-N0.4391x + 0.20810.50GNLSFSQI-A0.3088x + 0.29330.24GNLSFSQI-A0.3612x + 0.21660.44SQI-W0.3884x + 0.1510.39SQI-W0.1992x + 0.1850.11SQI-N0.3341x + 0.19810.24SQI-N0.3517x + 0.15590.42The MDS weighted additive soil quality index grades with the LLSF scoring technique represented more soil quality variation withR2 = 0.90 than the other soil qualities indices (Figure 2) (other spatial distributions of the other soil quality methods’ grade scoring techniques are not shown).Figure 2
Spatial distribution of MDS weighted additive soil quality index with the LLSF scoring technique.The results of sensitivity analyses and linear relationships between the indicator datasets with different scoring functions and index models point out that the MDS indicator dataset with the weighted additive quality model explains the soil quality indicators with the different quality models other than TDS and EDS indicator datasets and could assist in monitoring the temporal changes in soil quality of the region due to the increased agricultural development. These findings agreed with the findings of the other studies in arid and semiarid regions [36, 110, 121] where the MDS dataset with the weighted additive soil quality model is better for evaluating soil quality as the weighted additive model represents the indicators of the soil quality rather than the additive and Nemoro quality models [96, 108] and factor analyses provide the most representing indicators and the weights that differentiate the relative importance of soil quality indicators for evaluating soil quality [63, 107].


## 3.1. Soil and Terrain Indicator Characteristics

The soil and terrain indicators related to soil quality were subjected to descriptive analyses (Table2). The quantile range outlier’s statistic indicated no outliers existed in soil and terrain indicator data. The elevation of the sampling sites ranged between −19 and 197 m above mean sea level. The Rlf, sand, Slp, CaCO3, and clay indicators showed high data dispersion indicating that the spatial variations of these soil indicators are large, while the other soil and terrain indicators are of weak spatial variation as they have low data dispersion.Table 2
Descriptive statistics of soil and terrain indicators.
IndicatorMeanSEMinQ1MedianQ3MaxSand81.090.3965.0076.0080.0085.0099.87Silt6.160.170.095.007.5410.0017.79Clay8.510.280.048.0012.0015.0026.00pH7.380.027.007.167.247.508.20EC2.630.100.052.022.963.8613.14OM0.220.010.010.140.210.371.50ESP6.430.161.004.516.929.9813.46CEC2.130.100.021.462.323.1611.71CaCO31.950.130.041.442.634.2017.22FC10.950.232.899.5412.3915.0224.63AWHC5.400.062.864.895.606.358.61BD1.590.001.411.591.611.621.63PD1.690.001.421.671.721.751.84Slp5.900.220.544.096.329.8924.11Rlf14.940.283.0012.0015.0019.0044.00TWI7.790.095.526.827.248.3311.70Kaiser–Meyer–Olkin measure of sampling adequacyMSA0.74Bartlett’s test of sphericityApprox. chi-square16428.05Df120Sig.0.00Mean is represented by geometric mean, SE is the standard error, Min is the minimum value, Max is the maximum value,Q1 and Q3 are the lower and upper quartiles, respectively.The correlation analysis of soil physical and chemical indicators and terrain indicators revealed that there is a significant correlation in 35 of 136 soil and terrain indicator correlation pairs at a significant level ofp<0.05 (Table 3). The highest positive correlations were between soil clay content against soil field capacity, available water holding capacity, and soil packing density where the correlation coefficient is greater than 0.90. Also, the soil field capacity indicated a high positive correlation coefficient (r > 0.95) with available water holding capacity and soil packing density (r > 0.94). The highest negative correlation coefficients were obtained between sand and clay (r > 0.93), soil field capacity (r > 0.96), available water holding capacity (r > 0.99), and soil packing density (r > 0.93). It is also noted that soil organic matter indicator (OM) has a moderate positive relationship with CEC (0.45) and a moderate negative relationship with BD (−0.54) while it has no significant relationship with the soil salinity indicator (EC).Table 3
Correlation matrix for measured soil and terrain attributes (n = 397).
SandSiltClaypHECOMESPCECCaCO3FCAWHCBDPDSlpRlfSand1.00Silt−0.791.00Clay−0.930.511.00pH0.70−0.44−0.721.00EC−0.180.130.18−0.271.00OM0.08−0.09−0.060.31−0.031.00ESP−0.050.000.06−0.210.35−0.201.00CEC−0.040.030.040.270.020.45−0.081.00CaCO30.07−0.05−0.060.11−0.040.15−0.070.491.00FC−0.960.590.99−0.690.180.010.030.07−0.061.00AWHC−0.990.800.91−0.640.180.020.000.09−0.060.951.00BD−0.830.590.81−0.780.16−0.540.17−0.16−0.100.780.741.00PD−0.930.570.96−0.790.18−0.310.12−0.06−0.090.940.870.941.00Slp0.03−0.04−0.020.040.00−0.01−0.06−0.03−0.01−0.02−0.03−0.03−0.031.00Rlf0.04−0.01−0.040.10−0.020.08−0.110.080.06−0.03−0.02−0.09−0.060.511.00TWI0.06−0.13−0.010.000.06−0.130.09−0.14−0.12−0.03−0.080.040.02−0.14−0.04EC, electrical conductivity (dS m1); SOM, soil organic matter (%); ESP, exchangeable sodium percent (%); CEC, cation exchange capacity (cmol kg1); CaCO3, calcium carbonates content (%); FC, field capacity water content (%); AWHC, available water holding capacity (%); BD, soil bulk density; PD, soil packing density; Slp, slope (%); Rlf, relief (m), TWI, topographic wetness index.Bartlett’s test of sphericity revealed that the soil and terrain indicator correlation coefficients are significantly different than the identity correlation matrix (<0.05) and the indicators are not perfectly uncorrelated. The Kaiser–Meyer–Olkin Measure of Sampling Adequacy (MSA) revealed that the MSA value is above 0.6 and acceptable. This result of MSA indicates that there is a proportion of variance in soil and terrain indicators that might be caused by underlying factors. The Kaiser–Meyer–Olkin and Bartlett’s test show that there is a high scope for data reduction with factor analyzing where the pattern of data structure in indicators can be interpreted.


## 3.2. Minimum Dataset Indicators

The analyses of the relationship between soil and terrain indicators revealed that there are noted redundancy and collinearity among indicators and there is a need to establish a minimum dataset for the process of soil quality assessment by eliminating some indicators performing factor analyses.The results of factor analyses on standardized soil and terrain indicators values revealed that the first five factors accounted for 78.76% of indicator variability with eigenvalue greater than 1 (Table4). Thus, the first five factor analyses were retained and the other factors with eigenvalue less than 1 were excluded as the indicators may explain more variance than those excluded factors [117].Table 4
Eigenvalues, proportions and cumulative variance explained by factor analysis using a correlation matrix of standardized soil and terrain indicator data.
FactorEigenvalueDifferenceProportionCumulative16.726.5540.9440.9422.091.6510.3451.2831.521.539.5860.8641.281.459.0869.9451.001.418.8278.76The retained factors explain more than 97% of the variance for sand and AWHC in indicators (Table5); more than 91% for PD, FC, BD, and clay indicators; more than 71% for OM, CaCO3, Slp, CEC, pH, Rlf, and EC indicators; and more than 55% for ESP and silt indicators. The TWI indicator was the least important indicator due to its lowest communality estimate (less than 28% variance) as the higher communality estimates procure more preference upon the lower ones.Table 5
Proportion of variance using varimax rotation and communality estimates for soil and terrain indicators for each of the retained factors.
IndicatorFactor 1Factor 2Factor 3Factor 4Factor 5Communality estimatesSand−0.99−0.020.01−0.01−0.030.99Silt0.740.020.000.05−0.050.55Clay0.950.01−0.01−0.010.070.91pH−0.740.350.040.10−0.220.74EC0.170.130.05−0.070.810.71OM−0.090.920.030.09−0.080.87ESP0.00−0.23−0.09−0.020.790.68CEC0.030.55−0.030.670.070.75CaCO3-0.080.03−0.050.880.040.78FC0.970.09−0.01-0.010.050.96AWHC0.980.130.000.020.000.97BD0.85−0.47-0.040.000.100.95PD0.95−0.23−0.02−0.010.090.96Slp−0.02−0.070.880.02−0.020.78Rlf−0.030.090.850.06−0.020.73TWI−0.070.01−0.17−0.450.210.28The representative soil and terrain indicators for soil quality were chosen based on their relationship with the retained factors as the eigenvalues were the criteria magnitude for the interpretation (Table5). In the first factor, eight quality indicator factor loads are more than 0.7 which are sand, silt, clay, pH, FC, AWHC, BD, and PD. In the second factor, only the OM indicator is of high weight variable accounting for 0.92 of variation. In the third factor, slp and rlf indicators’ factor loads are of more than 0.85. CEC and CaCO3 indicators are more than 0.67 in the fourth factor. The fifth factor contains only EC and ESP indicator factor loads of more than 0.79. The soil and terrain quality indicators were grouped according to that criterion into five groups (Table 6).Table 6
Included soil and terrain indicators in the minimum dataset from representative indicators.
GroupIndicatorVector normCoefficient of determinationStandard normal transformationIndicator valueIncludedNDVIElevNormNDVIElev1Sand2.570.180.141.000.801.002.80No1Silt1.910.100.100.740.430.741.91No1Clay2.470.170.110.960.740.962.66No1pH2.010.020.010.780.110.781.67No1FC2.520.200.140.980.900.982.86No1AWHC2.540.230.170.991.000.992.97Yes1BD2.300.020.020.890.090.891.87No1PD2.480.090.070.960.390.962.32No2OM1.360.090.020.530.380.531.43Yes3Slp1.090.000.010.420.020.420.86Yes3Rlf1.060.000.010.410.010.410.83No4CEC1.100.000.010.430.020.430.87Yes4CaCO31.020.000.000.400.010.400.80No5EC0.950.010.000.370.030.370.77Yes5ESP0.860.010.060.330.020.330.69NoThe first factor explained about 41% of the total variance with high positive loadings from AWHC, FC, clay, PD, and BD indicators (>0.85) and high negative loading from the sand indicator (−0.99). The moderate positive loadings were from silt (0.74) and negative loading from pH (−0.74) indicators. The significant correlation coefficients among sand, silt, clay, pH, FC, AWHC, BD, and PD indicators influenced the factors’ loading values. The first factor expresses the influence of soil texture on hydraulic soil properties as it clarifies the soil moisture and water flow properties.The second factor explained about 10% of the total variance with high positive loading from the OM indicator (0.92) and moderate positive loading from CEC (0.55). The second factor had also moderate negative loading from the BD indicator (−0.47) and moderate positive loading from the pH indicator (0.35). The second factor expresses the effect of soil organic matter content on bulk density and soil cation exchange capacity.The third factor explained about 9.6% of the total variance with high positive loading from Slp (0.88) and Rlf (0.85) indicators. This factor is a proxy of surface as it clearly shows the effect of terrain on soil properties as they significantly correlated.The fourth factor explained about 8% of the total variance with high positive loading from the CaCO3 indicator (0.88) and moderate positive loading from CEC (0.67) indicators. Also, the fourth factor has negative moderate loading from the TWI indicator (−0.45) resulting from significant correlations among CaCO3, CEC, and TWI. The fourth factor expresses the interaction between CaCO3, CEC, and soil wetness as where wetness occurs, the calcium salts readily dissolved and interfere with soil exchangeable cations [118, 119] and the presence of calcium carbonates is associated with cascading changes in soil biogeochemistry [120].The fifth factor explained about 8.8% of the total variance with high positive loading from the EC indicator (0.81) and ESP (0.79). This factor represents salinity characteristics.The Norm analyses (Table6) revealed that, from eight quality indicators in the first group, only the AWHC indicator is with the highest Norm value. The OM indicator is the only indicator in the second group, while Slp, CEC, and EC indicators are the highest in the third, fourth, and fifth group, respectively. According to these results, the minimum dataset indicators include five soil and terrain quality indicators which are AWHC, OM, Slp, CEC, and EC indicators.


## 3.3. Soil Quality Assessment

The soil quality indices were calculated and ranked for each method. The classification criteria for the different calculated soil quality indices based on the different scoring techniques show variation between quality methods (A, W, and N), but there are slight similarities between each quality scoring technique for the different quality calculation methods. The threshold grade values of the scoring methods for each soil quality method were nearly the same threshold range values although the threshold grades values for the different scoring methods were not consistent.The cross validation for comparing the different selected interpolation ordinary kriging models identified the exponential predictive model as the best accurate model for representing the spatial distribution of soil quality indices other than the other predictive models (data not shown). The distribution of soil quality indices (Table7) shows a degree of pattern similarity between soil quality indices across the different quality grouping types with indicators scoring functions.Table 7
Distribution of index area grades across different soil quality index methods.
Indicator datasetsSQI modelStatisticSoil quality gradeI (very high)II (high)III (moderate)IV (low)V (very low)TDSSQI-AMin0.0013.8129.7312.340.12Max3.5231.0044.6935.823.53Mean1.7723.7039.6923.611.52SQI-WMin2.8018.0037.773.310.00Max5.7542.1949.7218.870.44Mean4.3135.1342.008.770.09SQI-NMin0.0016.8936.2612.110.11Max2.8425.3251.3728.081.52Mean0.8422.0246.6419.880.91MDSSQI-AMin0.0013.9319.480.170.00Max6.0863.0046.8933.350.17Mean2.4445.7433.828.260.03SQI-WMin0.0112.0429.750.400.00Max3.1258.4646.6236.550.08Mean1.0041.2238.099.960.02SQI-NMin0.001.9418.280.000.00Max3.9458.3247.2767.422.65Mean1.5337.4034.9915.830.55EDSSQI-AMin0.0013.2031.6813.110.00Max5.8828.8446.5237.504.04Mean1.5721.7840.5525.101.30SQI-WMin0.0012.5437.3812.640.00Max6.7828.3347.8128.816.08Mean1.7820.3742.9722.612.56SQI-NMin0.009.1438.7816.330.00Max5.6731.1553.0936.113.12Mean1.4322.6343.1222.410.70Soil quality grade III (moderate) is mainly the largest area percentage proportion in most soil quality indexing methods and scoring functions across the three soil quality datasets groups (TDS, MDS, and EDS) followed by grade II (High), grade IV (low), grade I (very high), and grade V (very low), respectively. A little part of the study area is for grades I (very high) and V (very low) with area percentage coverage. Soil quality grade III (moderate) on average accounts for 40% of the study area, and soil quality grade II (high) accounts for 30% of the study area, while the rest of soil quality grades (I, IV, and V) accounts for 20% of the study area. The disparities between the spatial distribution of soil qualities across indexing methods and scoring functions are of no large distribution disparity.


## 3.4. Comparison of Indices

Sensitivity analyses for soil quality indices showed that the soil quality grouping type of the indicator dataset, quality indexing method, and indicators scoring function are impacting the soil quality values in the study area (Table8). The LLSF linear scoring function is larger sensitive to the variation in soil quality compared with other sensitivities for all indicator grouping types, soil quality methods, and scoring functions specially with the minimum dataset grouping type while SNLSF nonlinear scoring function is the less sensitive to soil quality variation.Table 8
Sensitivity analyses between indicator datasets, models, and scoring functions.
Indicator datasetsSQI modelLSFNLSFLLSFHLSFGLSFSNLSFGNLSFTDSSQI-A2.481.771.621.392.14SQI-W3.683.191.851.532.68SQI-N2.521.891.821.692.14MDSSQI-A4.663.133.312.522.63SQI-W4.874.433.492.543.11SQI-N4.443.433.233.183.52EDSSQI-A3.612.582.212.023.38SQI-W3.943.722.462.422.83SQI-N3.383.063.282.193.09Sensitivity analysis showed that both indexing and dataset selection methods influenced the SQI values for the various soil quality values. Thus, the analyses emphasize that the indicator dataset grouping type, the indexing method, and scoring function could assist in the evaluation of the soil condition status in the study area. The Linear Scoring Functions (LSFs) disclosed more soil quality variations further than Nonlinear Scoring Functions (NLSFs). The resulted sensitivity analyses agree with the results in [89] where the Linear Scoring functions (LSFs) are more favorable than Nonlinear Scoring Functions (NLSFs) in evaluating soil quality indices. It is in contrast with the work in [15, 88] where it was reported that Nonlinear Scoring Functions (NLSFs) are favorable in studies with one linear and one nonlinear scoring functions, and their sensitivity was not evaluated [89].The linear relationships between the indicator datasets with different scoring functions and indices models (Table9) showed a significant correlation (P<0.05). The values of the linear relationships between MDS and TDS datasets are higher than the values of the linear relationships between EDS and TDS datasets across linear and nonlinear scoring functions and index models. The values of the linear relationships between TDS and MDS datasets for the weighted additive quality index are higher than the values of the linear relationship for additive and Nemoro indexing methods for all the linear and nonlinear scoring functions. In contrast, the linear relationships between the TDS and EDS datasets for the weighted additive quality index are mostly lower than the linear relationship values for additive and Nemoro indexing methods except for LLSF and HLSF linear scoring functions. The additive and weighted additive quality models’ values are quite similar for the scoring function type between the datasets.Table 9
The linear relationships between indicator datasets with different scoring functions.
SQIScoring methodSQI modelRelationshipR2SQIScoring methodSQI modelRelationshipR2TDS-MDSLLSFSQI-A0.4586x + 0.31320.54TDS-EDSLLSFSQI-A0.4634x + 0.23780.30SQI-W0.895x + 0.01080.90SQI-W1.0574x − 0.18690.64SQI-N0.5027x + 0.21110.51SQI-N0.4161x + 0.18650.23HLSFSQI-A0.3592x + 0.28880.46HLSFSQI-A0.4433x + 0.20450.43SQI-W0.61x + 0.06610.62SQI-W0.6876x − 0.01290.50SQI-N0.2621x + 0.19090.49SQI-N0.5026x + 0.13930.42GLSFSQI-A0.3725x + 0.24090.42GLSFSQI-A0.4703x + 0.150.39SQI-W0.4989x + 0.06650.51SQI-W0.5191x + 0.01430.31SQI-N0.4239x + 0.1610.41SQI-N0.5284x + 0.09970.39SNLSFSQI-A0.3487x + 0.3210.42SNLSFSQI-A0.3989x + 0.29960.50SQI-W0.4977x + 0.14890.48SQI-W0.4134x + 0.19260.32SQI-N0.3606x + 0.23190.39SQI-N0.4391x + 0.20810.50GNLSFSQI-A0.3088x + 0.29330.24GNLSFSQI-A0.3612x + 0.21660.44SQI-W0.3884x + 0.1510.39SQI-W0.1992x + 0.1850.11SQI-N0.3341x + 0.19810.24SQI-N0.3517x + 0.15590.42The MDS weighted additive soil quality index grades with the LLSF scoring technique represented more soil quality variation withR2 = 0.90 than the other soil qualities indices (Figure 2) (other spatial distributions of the other soil quality methods’ grade scoring techniques are not shown).Figure 2
Spatial distribution of MDS weighted additive soil quality index with the LLSF scoring technique.The results of sensitivity analyses and linear relationships between the indicator datasets with different scoring functions and index models point out that the MDS indicator dataset with the weighted additive quality model explains the soil quality indicators with the different quality models other than TDS and EDS indicator datasets and could assist in monitoring the temporal changes in soil quality of the region due to the increased agricultural development. These findings agreed with the findings of the other studies in arid and semiarid regions [36, 110, 121] where the MDS dataset with the weighted additive soil quality model is better for evaluating soil quality as the weighted additive model represents the indicators of the soil quality rather than the additive and Nemoro quality models [96, 108] and factor analyses provide the most representing indicators and the weights that differentiate the relative importance of soil quality indicators for evaluating soil quality [63, 107].


## 4. Conclusions

The objectives of this study were to (1) assess the soil quality using three types of indicator datasets (TDS, MDS, and EDS), three types of linear scoring functions, two types of nonlinear scoring functions, and three soil quality index models (additive, weighted additive, and Nemoro quality indices) and (2) find the most suitable indicator method and soil quality index model for the studied region with sensitivity and linear relationship statistical analyses.The findings of the study indicate the following:(i)
Factor analyses are an efficient tool to reduce the number of quality indicators for soil quality assessment, thus reducing the time and cost of sampling and analyses; also, it may improve soil quality assessment by increasing the sampling density. Exterior environmental variables are commensurate directly to soil quality. Factor analyses and exterior environmental variables identified five soil and terrain indicators, AWHC, OM, Slp, CEC, and EC indicators, to be included in the MDS dataset.(ii)
The linear scoring functions are more favorable than nonlinear scoring functions in reflecting the soil system functions. Although there was a variation between soil quality models (A, W, and N) based on the different scoring functions, for each scoring function, there were slight similarities across the quality models.(iii)
Grade III (moderate) soil quality grade is mainly the largest area percentage proportion across the three soil quality dataset groups (TDS, MDS, and EDS). The distribution disparity between resulted soil qualities is low with a degree of pattern similarity.(iv)
The MDS indicator dataset derived from factor analyses with the weighted additive quality model and LLSF scoring function was the most sensitive to assess soil quality as it is suitable to represent the TDS and explain the soil quality indicators.Based on the study findings, the researchers conclude that the selection of a minimum dataset of quality indicators in addition to inexpensive ancillary data is more beneficial, especially using remote sensing data and digital elevation models, for representing a sufficient soil quality index than using a more complex dataset. Moreover, the MDS approach is more advantageous in Egypt as a developing country where the requirements for soil quality assessment must be unpretentious and inexpensive to be adopted. The studied area is under reclamation and cultivated during the last two decades as a directive by the country to increase the agricultural areas and meet the increasing food requirements. Land degradation processes may be triggered as long-term cultivation continued coupled with improper land utilization. According to that, the results of the study could be significantly used as a starting strategy for future research to monitor the temporal changes in soil quality and utilization planning.



---

*Source: 1434692-2021-09-20.xml*