Soil Quality Assessment Based on a Minimum Data Set: A Case Study of a County in the Typical River Delta Wetlands

Abstract

Purpose: The Yellow River delta boasts rich land resources but lacks fresh water and exhibits poor natural conditions. To rationally develop and utilize the land resources therein, it is necessary to evaluate the soil quality. Methods: Adopting specific screening conditions, principal component analysis (PCA) was used to construct a minimum data set (MDS) from 10 soil indicators. Then, a complete soil quality evaluation index system of the Yellow River delta was developed. The soil quality comprehensive index (SQI) method was used to assess the soil quality in the Kenli District, and the soil quality grades and spatial distribution were analyzed. Results: (1) The average SQI of the Kenli District is 0.523, and the best soil quality is concentrated near the Yellow River, especially in Huanghekou town. (2) The normalized difference vegetation index was positively correlate with SQI, whereas Dr (nearest distance between the sampling site and Yellow River) and Ds (nearest distance between the sampling site and Bohai Sea) were negatively correlated with SQI. Elev (sampling site elevation) was not correlated with SQI. (3) The SQI of agricultural planting is greater than that of the natural land type and significantly greater than that of nudation. The main factors limiting farmland soil quality are SK (water-soluble potassium) and pH, whereas the factor limiting the natural land type are the soil nutrient indicators. Conclusions: To improve soil quality and develop and utilize land resources, the towns should adopt systematic land development/utilization methods based on local conditions. These results have important guiding significance and practical value for the more objective and accurate evaluation of soil quality in coastal areas and the development and utilization of land resources.

1. Introduction

Soil, air and water, are integral components of our environment [1,2,3]. Soil can support the growth of plants and animals, maintain and improve the quality of water and air, and provide adequate shelter for humans in both natural and unnatural ecosystems. Some scholars define soil’s ability to achieve these functions as soil quality [4]. Understanding soil quality is important to improving sustainable land use management [5], providing early warning signals of adverse trends in soil quality change, identifying problematic areas of soil quality [6], and providing a valuable basis for the subsequent rational use and improvement of soil [7]. However, these understanding can only be obtained through reliable and accurate soil quality evaluation methods. Soil quality evaluation refers to the monitoring and evaluation of soil properties, soil functions and soil conditions.

Accurate evaluation of soil quality is difficult given the heterogeneity and variabilities in the physical, chemical and biological characteristics of soil in different regions [8,9,10], and the use of pesticides and fertilizers and different soil management measures further complicate the evaluation [11]. In addition, the soil quality evaluation is often confused with soil fertility evaluation. Thus, in the process of a soil quality assessment on soil fertility status [12], the limiting factors of soil quality vary according to the type of land use, the geographical location of the ecosystem and the soil parent material [13,14]. At present, there are no recognized norms or standards used to establish soil quality evaluation index systems or soil quality classifications [15,16]. Based on the variability and complexity of soil properties, a variety of soil quality evaluation methods have been developed, such as the soil quality index (SQI) [17,18], fuzzy association rules [19], the dynamic soil quality model [20], the soil management evaluation framework [21], and the soil quality card and detection kit [22]. Among these methods, the SQI is widely used because it can be applied to flexibly quantify different soil types and can flexibly integrate the physical, chemical and biological properties of soil [15]. The quantitative calculation method of soil quality index adopts three steps: indicator selection, indicators scoring and comprehensive score into one index [14,23,24].

By establishing the minimum data set (MDS), we can easily select the most appropriate indicator to evaluate soil quality, and this indicator has been used widely [4,25,26]. The selection of MDS indictors can reduce the redundancy of soil data. In addition, the weight of selected indicators can be calculated while establishing the MDS, thus reducing the subjective influence [27]. MDS can be screened by linear and multiple regression factor analysis, discriminant analysis and score function. Because a factor analysis can reduce redundant information in the original data set, it is widely used in the identification of the MDS [28]. For example, Peng Li used a one-way analysis of variance (ANOVA) and factor analysis to select four of 26 soil indicators and built an MDS to evaluate soil quality in different wheat-producing areas of China [4]. Research on soil quality assessments based on MDS have been effectively applied and promoted in different land use types, such as offshore areas [28,29], agricultural areas [8,15] and grasslands [27]. Considering the external environment and human management influence on soil quality, some scholars have attempted to add macro soil environmental factors and land use status into the indicator selection criteria to improve the establishment of the MDS and have achieved good results, however, this method has not been widely used [30,31].

The Yellow River delta was mainly formed by sediment deposition from the Yellow River [32]. A large number of nutrients are carried by the sediment in the upper reaches of the Yellow River. However, these nutrients are simultaneously affected by seawater intrusion, and soil salinization is becoming increasingly serious [33]. In addition, in recent decades, there have been frequent human activities such as oil extraction, cultivation and urban construction in the Yellow River delta [32,34]. Moreover, water inflow from the upper reaches of the Yellow River decreases every year, and coastal wetlands have significantly degraded [34]. Influenced by natural conditions and human activities, natural ecosystems have been degraded by artificial isolation [32]. Overall, the soil quality has deteriorated, but in some areas, soil quality has been improved by tillage, which results in the great variation in soil quality in the Yellow River delta [31]. These problems in the Yellow River delta hinterland, the Yellow River estuary Kenli District, are particularly obvious. In recent years, studies of soil quality in the Yellow River delta have been relatively limited. In addition, some problems were noted with previous studies. The Yellow River and the ocean are important regional features of the Yellow River delta; however, previous studies on soil quality in this region have not considered the impact of external factors, such as the Yellow River and coastline, on soil quality [35]. Moreover, relevant studies are not comprehensive or scientific in nature.

The objective of this study was to evaluate the soil quality of the typical coastal area of the Yellow River delta-Kenli District based on the MDS method. The vegetation normalization index, the distance from sampling points to the Yellow River and the coastline and elevations of sampling points were considered the influential factors in the indicator screening process. This study attempted to find out the limiting factors of the soil quality characteristics of different land use types and towns, as well as the correlations between the external environment factors and soil quality, some of which were firstly studied in the Yellow River delta. The results of this study have important guiding significance and practical value for objectively and accurately evaluating the soil quality of reclamation areas in coastal areas.

2. Materials and Methods

2.1. Study Area

The Yellow River delta wetland is the largest, most complete and youngest estuarine wetland ecosystem preserved in the warm temperate zone of the world. The delta is formed by the interaction of the Yellow River alluvium and marine deposition, with loose lithology and low and flat terrain. The region has a subhumid continental monsoon climate with four distinct seasons and rainy summers. The average annual temperature is 12.9 °C, and the average annual precipitation is 560 mm, 70% of which occurs from July to September [33]. The average annual evaporation rate is 1900–2400 mm, with an evaporation precipitation ratio greater than 3.5 [32]. Affected by tides and Yellow River runoff [36], the average groundwater level in this area is 1.1 m [32], and the average groundwater salinity is 30.1 g/L [37].

The Kenli District in Dongying city of Shandong Province was selected as the research area. The Kenli District is located at the estuary of the Yellow River in the hinterland of the Yellow River delta. The regional location is 37°24′–38°10′ N, 118°15′–119°19′ E. Due to its low elevation (generally less than 10 m except for roads and levees) and proximity to the ocean, the hydrological characteristics of the Kenli District are affected by the interactions between freshwater and seawater, as well as between groundwater and surface water [34]. Due to the large evaporation precipitation ratio, high groundwater level and salinity, soil salinization in the study area is a more serious problem. Repeated breaches that have been diverted from the Yellow River have resulted in the regional formation of hillocks, slopes, depressions and landforms. In addition, due to human activities, such as the artificial diversion of the Yellow River, Yellow River levee and coastal levee building, oil exploitation, urban construction, and farming techniques allowed on the delta, great changes have occurred in these micro landforms. These changes determine surface runoff and groundwater movement, forming the water and salt gathering area with the depression as the center and a “post drought, hollow waterlogging and the second slope basification” phenomenon. Due to its unique geographical location, in a normal year, the Yellow River can produce approximately 20 km² of new land by carrying sand. The Yellow River sediment brings a large amount of nutrients from the upstream to the Kenli District. The main natural vegetation in this area is reeds (Phragmites australis), Tamarix chinensis, and Spartina alterniflora. Artificial vegetation mainly includes cotton, corn, rice, wheat, legume and white wax.

2.2. Data Collection

Before the on-site sampling was conducted, a total of 149 sampling points were set up evenly in the Kenli District and nearby areas at the shortest distance of 2 km by using Google Earth software and the visual interpretation method. Then, the WGS84 coordinates of the designed sampling points were imported into a handheld GPS navigator. Field sampling was conducted according to the predetermined location of sampling points and the field conditions. Soil drills with a diameter of 5 cm were used to collect the surface soil at 0–20 cm in 3 places within a range of 3 m from each point. After mixing, the samples were sealed in bags and numbered. A handheld GPS navigator was used to record the exact longitudes and latitudes of the sampling points, and the land use types, vegetation, hydrological characteristics and topography around the sampling points were also recorded. Finally, a total of 127 sampling points fall within the scope of the Kenli District (point position distribution is shown in Figure 1), and a total of 9 planting cover types were sampled in 7 townships (Table 1). Because most of the soil is flooded by the seawater, there are no sampling points in the intertidal zone, and the study area does not include the intertidal zone. To ensure the reliability of the statistical analysis, a comparative analysis is not conducted for vegetation types with fewer than 6 sampling sites, and a statistical comparative analysis is not conducted for towns with less than 8 sampling sites.

2.2.1. Soil Sample Analysis

The collected soil samples were naturally air-dried in the laboratory, and the roots and other sundries were removed. Then, the soil samples were ground and sifted through a 1-mm sieve. The ground soil samples were accurately weighed to 6.0000 g, and placed in a centrifuge tube. Then, 30 mL of deionized water was added. After oscillating on the oscillator for 5 min, the samples were centrifuged at 10,000 r/min for 5 min. Then, the supernatant was extracted to prepare the soil extract with a soil to water ratio of 1:5. Partial soil chemical parameters were determined using the method described in the literature [38]: The pH of the soil extract was determined by a PHSJ-4A pH meter, and the conductivity of the extract (EC) was determined using a DDSJ-308A conductivity meter. The soil organic matter (OM) was determined using the potassium dichromate external heating method. The soil extract obtained by centrifugation was diluted by 2 mL 20 times and then sent to the testing and analysis center of the Yantai Institute of the Coastal Zone, Chinese Academy of Sciences for ion concentration determination. The concentrations of NH₄⁺, NO₃⁻, Cl⁻, Na⁺, SO₄²⁻, K⁺, Ca²⁺, and Mg²⁺ in solution were determined by a Dionex ICS3000 ion chromatograph. The sum of the mass of NH₄⁺ and NO₃⁻ was used to represent the soil available nitrogen (AN), and K⁺ was the soil water-soluble potassium (SK).

2.2.2. The Auxiliary Data

The auxiliary data are used as the influencing factor of the external environment in the indicator screening process [39,40], including the altitude of the sampling point (Elev, the data are obtained from Google Earth geodesic elevation, and the sampling interval is 7.57 m), Normalized Difference Vegetation Index (NDVI, data were extracted from Landsat 8 OLI at a sampling interval of 30 m), the nearest straight-line distance between the sampling point and the Yellow River (Dr, data were obtained by using ArcMap c after linear vectorization of the Yellow River), the nearest straight-line distance between the sampling point and the coastline (Ds, data were obtained using ArcMap nearest neighbor analysis after linear vectorization of the coastline) [41]. All data sets are projected into a universal transverse Mercator Cartesian projection coordinate system based on WGS84 (the EPSG code is 32650).

2.3. Methods of Soil Quality Assessment

Soil quality index was calculated by four steps (Figure 2). Firstly, correlation analysis is carried out to determine whether the soil indicators were redundant. Secondly, principal component analysis and external environmental indicators are combined to screen the minimum data set (MDS) and obtain their weights. Thirdly, standard scoring functions (SSF) were used to score the MDS. Finally, the soil quality index is calculated according to the indicator score and its weight. The specific processes are described as follows.

2.3.1. Determine the MDS

Soil quality assessment is based on the selection of representative soil indicators. We used principal component analysis (PCA) among soil indicators and linear regression coefficients between them and external environmental factors to determine the appropriate MDS for soil quality assessment [31].

First, the grouping of evaluation indicators was determined. A Pearson correlation analysis was performed on soil indicators to analyze whether collinearity existed between indicators and determine whether it was necessary to screen out all the indicators from preselected indicators. If there is redundancy in the soil indicators, a PCA is used to group the indicators, with the component of the eigenvalue ≥1 as the main component. Indicators in the same principal component with a factor load of ≥0.5 were placed into one group. If the factor loads of a certain indicator in different principal components are all ≥0.5, the indicator was placed into the group with the lowest correlation with other indicators. If the correlation coefficient between an indicator and other indicators in a group is less than 0.3, then the indicator is placed into a separate group. Here, the indicator grouping has been determined and completed.

Then, the MDS is determined by calculating the indicator values of each indicator. The indicator vector norm is calculated as noted in Formula (1). Considering that soil quality is determined by both internal properties and the external environment [42], the elevation of sampling points (Elev), NDVI, distance to the Yellow River (Dr) and distance to the coastline (Ds) were added as the screening criteria in this study. Through a regression analysis, the coefficient of determination (R²) between the factors influencing the external environment and each indicator was obtained, and R² reflected the degree of influence of the external environment on the soil indicators. Then, the vector norm of each indicator and the determination coefficient R² are summed by the standard linear transformation to obtain the indicator value. The indicators set to at least 90% of the maximum indicator value in each group are the preselected indicators. If the correlation coefficient of any two preselected indicators is greater than 0.5, then the indicator with the higher indicator value is selected for the MDS. Otherwise, all metrics are selected as MDS.

(1)$N_i=\sqrt{{\displaystyle \sum _{j=1}^k}{u_{ij}}^2{\lambda }_j}$

$N_i$ is the combined load of indicator i in all components with eigenvalues ≥1,$u_{ij}$ is the load of indicator i in principal component j, and ${\lambda }_j$ is the eigenvalue of the principal component j.

2.3.2. Weight Allocation

After the MDS indicator is determined, the weight of each MDS indicator is the ratio of the common factor variance in the PCA to the sum of the common factor variance in all MDS indicators.

2.3.3. Indicator Scoring

After determining the indicators of the MDS, each indicator was converted using a standard scoring function (SSF), which was normalized to a value between 0.1 and 1 according to the characteristics of the soil indicators [43]. Generally, the indicators’ scoring functions follow “more is better”, “less is better” or “optimum” scoring curve trends. Equations (2)–(4) were used for scoring of “more is better”, “optimum” and “less is better” curves, respectively:

(2)$SSF1:f\left(x\right)=\{\begin{array}{c}0.1, x<L\\ 0.1+\frac{0.9\times \left(x-L\right)}{\left(U-L\right)}, L\le x\le U\\ 1.0, x>u\end{array}$

(3)$SSF2:f\left(x\right)=\{\begin{array}{c}\begin{array}{c}0.1, x⟨L_1, x⟩U_2\\ 0.1+\frac{0.9\times \left(x-L_1\right)}{\left(U_1-L_1\right)}, L_1\le x\le U_1\\ 1.0,U_1< x<L_2\end{array}\\ 0.1+\frac{0.9\times \left(U_2-x\right)}{\left(U_2-L_2\right)}, L_2\le x\le U_2\end{array}$

(4)$SSF3:f\left(x\right)=\{\begin{array}{c}1, x<L\\ 0.1+\frac{0.9\times \left(U-x\right)}{\left(U-L\right)}, L\le x\le U\\ 0.1, x>U\end{array}$

where x is the measured value of the soil indicator, $f\left(x\right)$ is the score of the soil indicator, and U and L are the upper and lower critical values of the soil indicator, respectively.

2.3.4. SQI Calculation Method

Finally, SQI was calculated according to Formula (5) [44]:

(5)$\mathrm{SQI}={\sum }_1^n(W_i\times S_i)$

where $W_i$ is the weight of indicator i; $S_i$ is the score of indicator i, which was calculated according to SSF; and n is the indicator number of the MDS.

2.4. Statistical Analysis

SPSS 25.0 was used for the statistical analysis. The Kolmogorov-Smirnov test (K-S normal test) was used for data normality testing. The Pearson product difference correlation test (bilateral) was used for the correlation analysis. PCA analysis was used to group soil indicators. A one-way ANOVA was used to quantitatively analyze the differences in SQIs of different vegetation types and different townships. A linear regression is used to quantify the influence of external factors on soil indicators and soil quality. To test the interpolation accuracy, 16 methods were used in this study: simple to quartic inverse distance weight (IDW), global first degree polynomial (GP), simple to quartic local polynomials (LP), tension spline interpolation (SPT), regular spline interpolation (CRS), higher-order surface function interpolation (MTQ), inverse higher-order surface function interpolation (IMTQ), thin plate spline function interpolation, simple kriging (SK) and ordinary kriging interpolation (OK). The interpolation diagram was drawn by Golden Software Surfer 15.0 [45], and the other figures were generated using Python 3.7 Pandas, matplotlib, and seaborn.

3. Results

3.1. The Characteristics and Spatial Variation Description of Soil in the Kenli District

Ten soil characteristics and four external environmental factors related to soil quality were measured and analyzed (Table 2). Values of pH, OM and all four external environmental factors (NDVI, Elev, Ds, and Dr) passed the K-S normal test. From the spatial position of the sampling points, the nearest points to the Yellow River and the coastline are only 464 and 1040 m, respectively, and the farthest are 24.8 and 55.7 km, respectively. The elevations of the sampling points range from −0.94 to 12.37 m. The average OM was 10.15 g/kg, and the maximum OM was 21.77 g/kg. The average pH of the soil was 8.50, which was considered to be a weakly alkaline soil with a weak spatial variation. Other soil indicators (EC, AN, Cl⁻, Na⁺, SO₄²⁻, Ca²⁺, Mg²⁺, and SK) are lognormal, and the large standard error and high data dispersion indicate that the spatial variation in these soil indicators is large. Table 2 demonstrates that the maximum value of the indicator is much greater than the upper quantile, indicating a high degree of data dispersion.

3.2. MDS Selection and SQI Calculation

3.2.1. MDS Selection

The Pearson correlation coefficients of all soil character indicators (Figure 3), including EC, Cl⁻, Na⁺, SO₄²⁻, Ca²⁺, and Mg²⁺, were significantly correlated with each other, and the correlation coefficient was greater 0.7. The correlation coefficient among EC, Cl⁻ and Na⁺ was greater than 0.97, which indicates that Cl⁻ and Na⁺ are the most important ions affecting soil salinity and that seawater salinity is the main source of soil salinity. The pH was negatively correlated with all the indicators, and the correlation coefficients were the highest with Mg²⁺ and Ca²⁺, which were 0.63 and 0.73, respectively, indicating that it was feasible to add calcium and magnesium modifiers to reduce the pH of the saline soil. OM and AN are not significantly correlated with most soil indicators. OM is only significantly negatively correlated with pH and Cl, but the correlation coefficient is low (0.24 and 0.17, respectively), indicating that soil salinity has no significant relationship with OM. AN is only significantly correlated with SK, Ca²⁺ and pH, but the correlation coefficient is relatively low at between 0.1 and 0.25. We divided the clustering results into five categories: “SK”, “OM”, “AN”, “EC, Cl⁻, Na⁺, SO₄²⁻, Ca²⁺, Mg²⁺”, and “pH”. The classification results represent soil fertility properties (potassium, nitrogen, and organic matter, which are the first three categories), salinity characteristics (conductivity and salt ions, the fourth category) and acid-base properties (pH, the fifth category). The correlation analysis results reveal that EC, Cl⁻, Na⁺, SO₄²⁻, Ca²⁺, Mg²⁺ and SK are collinear with each other, and redundancy is noted in the interpretation of soil properties. Therefore, it is necessary to eliminate some indicators in the process of establishing the MDS.

Further, we used the PCA to select an MDS from 10 soil indicators [44]. Three principal components (PCs) that had an eigenvalue of >1, and the cumulative variance interpretation rate is 80.226%. In PC1, there are 7 indicator factor loads greater than 0.5, which are mainly related to soil salinity properties. In PC2, pH, AN and OM were the high weight variables, accounting for 26.68% of the variation. In PC3, AN and SK are high weighted variables that represent the nutrient status of soil. Then, all the indicators were divided into five groups according to the above grouping criteria.

The vector norm values (norm) of each indicator are calculated according to Formula (1) and the eigenvalues and factor loads are listed in Table 3. The Pearson correlation coefficient between the soil indicators and external environmental factors (Elev, NDVI, Dr, and Ds) was obtained by SPSS correlation analysis (Table 4). The analysis results demonstrate that EC, Na⁺, Cl⁻ and SK, the four ions that can represent soil salinity, are significantly correlated with the four external environmental factors: they are positively correlated with Elev, NDVI and Ds and negatively correlated with Dr. No significant correlation was noted among AN, OM and external environmental factors, indicating that the distance between the elevation of vegetation coverage and the coastline of the Yellow River did not affect the distributions of AN and OM. pH was significantly positively and negatively correlated with Ds and Dr, respectively, with a correlation coefficient of 0.24. However, this correlation was not related to NDVI or altitude, indicating that seawater would increase the soil alkalinity, whereas fresh water could improve the soil alkalinity [46]. Except for OM and SK, the correlation coefficients of Ds, Dr and other soil indicators are very similar, indicating that the influence of seawater and the Yellow River on most of the soil chemical properties is negative and equivalent. No significant correlation exists between the indicator and the external environment, the coefficient of determination is 0. Finally, according to the above MDS screening method, the MDS included the five indicators of Na, OM, pH, AN, and SK (Table 5).

3.2.2. Calculation of the SQI

After the MDS was determined, the soil quality was evaluated through three steps. First, the weight of the MDS indicator is calculated by the communality derived from the PCA (Table 6). The results show that Na and AN are more important than pH, OM and SK in the soil quality evaluation because they have a higher weight. Furthermore, according to the grading standards of the second soil investigation in China and the relevant literature [31], the upper and lower thresholds of soil MDS evaluation indicators are determined. According to the impact of indicators on soil quality, the standard linear conversion function with the lower limit type described in Formula (2) was used for the standardized scoring of AN, OM and SK, whereas the opposite, Na were scored using the standard linear conversion function described in Formula (4). In general, pH was scored using the one described in Formula (3), however, the soil in the study area was slightly alkaline with a pH greater than 7, so we scored it use Formula (4). Finally, the SQI is calculated by Formula (5).

3.3. Soil Quality Assessment in the Kenli District

3.3.1. The SQIs of Different Land Use Types

The SQI of the 127 sampling points ranged from 0.237 to 0.836, with an average value of 0.523, a coefficient of variation of 26.0%, and a moderate intensity variation. We used the soil quality classification standard in Wu’s paper [31] for the Yellow River delta to classify the soil quality grade of the study area. According to the SQI > 0.7, 0.6–0.7, 0.5–0.6, 0.4 to 0.5 and ≤0.3, the soil quality can be divided into levels I to VI, and I is the best level in terms of soil quality. The SQI of all sampling points was graded, and the results (Table 7) demonstrated that the proportion of grade IV soil was the highest, followed by grade III soil. The total proportion of these two grades was 52.7%, whereas the proportion of grade VI soil was the lowest at only 3.9%. The grade I point number was 2.6-fold that of grade IV.

As evident in Figure 4A, the mean SQI of agricultural land (cotton and corn) was significantly increased compared with that of nudation, and the mean SQI of natural land type (grassland, forest, and reed) was lower than that of agricultural land and higher than that of nudation, without significant differences. The spatial variability in the SQI of farmland type is significantly increased compared with that of the natural land type and bare land type under the influence of cultivation years and cultivation modes. Due to its good adaptability to the environment [47], reeds have a wide spatial distribution and a large variation coefficient in the SQI. This feature is also noted in the distribution of soil quality grade. The soil quality of the cotton sampling points span across all the soil quality grades. Corn simultaneously also has grade I soil and grade VI soil, as well as various soil quality grades across a large span. The soil quality grade of the natural type is mainly grade IV, and the soil quality of nudation is worse than grade IV (mainly grade V).

According to the indicator contribution value and contribution rate of each land type’s SQI (Figure 4B), in agricultural systems, the pH is high, SK is low, and the main limiting factors of farmland soil quality are SK and pH. Natural type and nudation generally have lower AN and OM values, which are the main factors limiting soil quality. In addition, high Na was also an important factor limiting the soil quality of grassland and pH of forest. Compared with natural vegetation types and agricultural systems, Na, which characterizes soil salinity, contributes minimally to SQI on bare land. The contribution rate of grassland pH to the soil SQI is significantly increased compared with that of other types.

3.3.2. Spatial Distribution Characteristics of Soil Quality

In Figure 5A, differences in the soil quality SQIs are noted between different towns in the Kenli District. Among them, the value for Huanghekou town is significantly greater than that for Yong’an, Dongji and Shengtuo. The value for Kenli Street is less than (not significantly) that for Huanghekou town but greater than (not significantly) that for the other three towns. No obvious spatial pattern is noted.

There are differences in the contribution of the MDS indicator to the SQI in each town (Figure 5B), reflecting the differences in soil ecology and soil management in different towns: Yong’an town is close to the coast and far from the Yellow River, and the amount of soil salt is very high. Kenli Street pH is high, which is the main limit on soil quality. No significant difference in the contribution of soil AN to soil quality was noted among the townships, but the value of the two towns in the west was lower than that on Kenli Street and the two towns in the east. The SK of Shengtuo town is significantly lower than that of the other townships, so this area should pay attention to the addition of potash fertilizer in its land use. Supplemented by the high OM in the alluvial sediment from the Yellow River, the contribution of OM to soil quality in Huanghekou town was the largest among all the townships and significantly greater than that in Dongji town and Yong’an town.

Then, we performed a geostatistical analysis of the SQI of the study area. The results of the spatial autocorrelation analysis show that Moran’s I index = 0.13, z-score = 1.65, and the probability of randomly generating this clustering pattern is less than 10%, which means that the SQI is spatially autocorrelated. Among all the models, the linear model can best explain the spatial variability in the SQI (R² = 0.525, RSS = 1.90 × 10⁻⁵); the nugget effect was 0.016, and the base station value was 0.020 (Figure 6).

By comparing 16 interpolation methods, IMTQ was identified the best method for spatial interpolation of the SQI (RMSE = 0.127), and the worst interpolation method was TPS (RMSE = 0.166). Figure 7 presents the IMTQ interpolation diagram of the spatial distribution of the SQI grade in the study area.

Figure 7 demonstrates that the spatial distribution of soil quality has a strong regularity. Generally, the closer the study area is to the Yellow River, the better the soil quality is. Grade I soil is mainly concentrated in the vicinity of the Yellow River, especially in Huanghekou town. The distribution of grade VI soil is relatively scattered. Due to the high degree of urbanization around Kenli town, the soil quality is relatively poor. Overall, the soil in the Kenli District is poor in the central and western regions, whereas the eastern part has the best soil quality. According to the statistical analysis, the proportion of soil quality grades I–VI in the study area was 4.5%, 13.2%, 33.0%, 42.2%, 6.1% and 1.0%, respectively. Grade III and grade IV account for the largest proportion, 75.2% of the total area, which is 22.5% greater than the statistical result of the sampling points. Grade I soil accounts for 4.5 times more area than grade VI soil, which is greater than the statistical result of the sampling points.

3.4. Analysis of External Environmental Influencing Factors of SQI

The ocean, Yellow River, microtopography and vegetation change play important roles in land generation, soil morphology, soil nutrients and salinity redistribution. In this study, the SQI was used as the dependent variable, and NDVI, Dr, Ds and Elev were used as independent variables for the linear regression analysis (Figure 8).

The correlation coefficients between the four selected external environmental factors and SQI were relatively small. However, SQI and NDVI exhibited a very significant positive correlation (Figure 8A), and SQI was somewhat significantly (0.05 < p < 0.1) negatively correlated with Dr and Ds (Figure 8B,C). However, no correlation was noted between Elev and SQI (Figure 8D). Although the elevation of the study area fluctuated, the overall slope was approximately 1/10,000, which was not sufficient to significantly affect the soil quality.

4. Discussion

4.1. Scientific Evaluation and Explanation of the MDS

The MDS method is widely accepted in the literature because it is carefully selected from a wide range of soil indicators [27,48]. In this paper, the MDS included Na, pH, AN, OM and SK, which were screened by combining external environmental factors and the PCA method. Indicators used in the soil quality assessment should be site-specific and plant-specific [49]. The MDS indicators obtained in this paper fully considered the comprehensive soil chemical indicator and its effect on crop ecological plant growth in the Kenli District of the Yellow River delta. Some of these indicators (pH, AN, OM, and K) have also been widely used in many studies [4,28,29,31]. Wu et al. selected the minimum data set including TN, pH, AK and SOM to evaluate the soil quality of the Yellow River delta using the principal component analysis method assisted by environmental factors [31]. Zhang et al. used PCA to select AN, SOM, MBN(Microbial biomass carbon)/TC, and available zinc from 26 soil indicators as the minimum data set to evaluate the soil quality of the wheat-maize cropping system in the region of the dryland double cropping, China [4].

Nitrogen is one of the most important soil nutrients in farmlands and is a limiting nutrient in most ecosystems [50], especially effective nitrogen [51]. Soil OM is an important component of many soil biophysical properties [52], and plays a key role in improving soil fertility and ecosystem productivity [53,54,55]. Potassium is one of three essential nutrients for plant growth [56], and plays an important role in plant growth and metabolism. In addition, potassium can alleviate stress in plants under drought and salinization conditions [57]. Therefore, soil available potassium should be used as an important indicator of comprehensive SQI evaluation. Soil pH is an important indicator of soil alkalization.

Na is rarely selected as an MDS indicator in soil quality evaluations. However, in this paper, Na was selected as the characterization indicator of soil salinity in the Yellow River delta rather than EC because the former could reflect the influence of the external environmental factors on soil quality better than the latter (Table 4). Soil salinity results in growth inhibition, developmental changes, metabolic adaptations, and ion sequestration or exclusion. Both osmotic and toxic effects of sodium ions are involved [58]. These two processes were thought to be temporally and spatially separated, with sodium ions in the soil rapidly reducing water availability and a slow accumulation of sodium in the shoot inhibiting photosynthesis [59].

The assessment result and indicator screening were reasonable, as the spatial distribution of soil quality was similar to the actual conditions that we surveyed in the field. However, there was a flaw in this study. Because of the restrictions of experimental conditions, we could not get and analyze the characteristics of biological and physical indicators and their roles played in soil quality, and this was what we would solve in later researches.

4.2. Rationality of the Soil Quality Evaluation in the Study Area

According to Wu et al. [31], garden plots had the highest mean soil quality (value is 0.61) and the land use type with the second highest soil quality was farmland with a mean value of 0.58, and farmland primarily had grade III soils. In our results, although corn fields have the best soil quality among all land types, the average SQI is just 0.60 (Table 7) and is classified as grade III, which are consistent with Wu’s. In a similar study, Zhang et al. applied the same methodology for evaluating the SQI of the Yellow River delta wetland, they reported an SQI less than 0.53 and soil salt content had the highest SQI contribution ratio (34.1%) among the indicators [29], which is exactly similar to those found in our study. Soil salt content was the limiting factor of soil quality due to its high contribution ratio to the SQI value (Figure 4A and Figure 5A), indicating that the reclamation area in the Yellow River delta has a high salinity and a relatively poor soil quality. Under the conditions of artificial fertilization and soil mechanical improvement in saline land, the soil quality of agricultural land (cotton and corn) was greater than that of the natural land type (grassland, forest, and reed), and the soil quality of the natural land type was greater than that of nudation. Meng et al.’s research in east coastal region of China has come to the same conclusion, changes in different land use types due to plants and application of fertilizers were characterized by promoted soil quality [60]. And Gong et al.’s research for the saline oasis soil of the Keriya River Basin also reached a conclusion that concentrations of soil nutrients increased and land improvement occurred after land use shifts from natural land uses to cultivation [61].

Over the past few decades, increasing anthropogenic activities, such as urbanized agricultural activities and oil exploration, have led to rapid changes in land use and ecosystem fragmentation in the Yellow River delta. Since the 1980s when China implemented the household contract responsibility system, the area of arable land operated by individual farmers has been very small [62]. The spatial variability in soil nutrients is much greater than that in large-scale agriculture in developed countries [62]. These are the important factors that cause differences in the SQIs between different towns and the soil at all grades dispersed throughout the Kenli District. Nevertheless, due to the fresh water irrigation of the Yellow River and the rich nutrition of the new land near the coastal zone [63,64], the high quality soil was mainly distributed in areas close to the Yellow River and the coastline. The analysis of the influencing factors also showed that the Yellow River and coastal zone could promote the soil quality. Moreover, Wu’s results showed that grade I was distributed only near the crossing of the Yellow River and the Diao River, which also supported our results [31]. However, the spatial variation and irregularity of the soil quality in this paper are significantly greater than those in Wu’s study because the sampling density in this paper is higher, and the intertidal zones without sample points are eliminated in the interpolation process, the interpolation analysis may be more elaborate.

4.3. Recommendations for Soil Management in the Yellow River Delta

In the process of transformation of original natural land use to farmland, land degradation usually occurs [65,66]. Nevertheless, this conversion of natural vegetation to farmland does not always lead to land degradation [67,68]. Agricultural planting and natural vegetation growth could significantly reduce the soil salinity and improve the soil quality of saline-alkali soils [60,61], and a large number of studies have demonstrated that vegetation and soil quality can support each other [69]. Thus, in the process of developing and utilizing the coastal saline soil of the Yellow River delta, we can reasonably develop wastelands or protect natural vegetation systems to prevent degradation and nudation. Figure 4B shows that the soil pH can be reduced and the ecological environment can be improved by planting grass for high pH soils that are not suitable for development and utilization. In agricultural systems, pH is high, SK is low, and potash should be continuously added to improve soil quality and ensure crop yield. Figure 5B shows that the limiting factors of soil quality were different in each township due to the influences of the environment and the geographical location. In view of the differences and limiting factors of the soil quality of different land use types in different towns, we should propose land development and utilization strategies. For example, the pH of Kenli Street is high, and this area is highly urbanized, and the results above show that grassland can reduce the soil pH, which is a limiting factor for the soil quality of forests. Therefore, in the process of soil greening in the Kenli Street area, we should consider planting lots of grass and planting trees carefully.

The formation and evolution of the land in the Yellow River delta is greatly influenced by the external environment [70]. The interaction between the Yellow River sediment and the Bohai Sea tide resulted in the formation of high-salt and high-nutrient land, and the plant growth stabilized the soil and improved the saline soil [34]. The results of this paper also proved that NDVI was directly proportional to soil quality. The Yellow River and Bohai Sea as the main material sources have the most obvious influence on the soil quality of the Yellow River delta. The soil near the Yellow River has high nutrient content, abundant fresh water resources and better soil quality. The soil in this area is good for agriculture. Close to the ocean, the soil is salty, but the tides bring nutrients. So, the area could be developed for salt-water farming or for salt-tolerant crops. These results are of great significance to soil development and utilization planning.

The study indicated that reasonable land use practices and scientific management patterns can improve and maintain the soil quality in the human-disturbed land-use systems in the Yellow River delta. However, other research has demonstrated that long-term cultivation can lead to soil degradation [71], which may explain the poor soil quality in the two western townships. Given the vulnerability of the Yellow River delta ecosystem, unreasonable land use can aggravate land degradation processes within fragile ecosystems. Therefore, sustainable land use and management are particularly important.

5. Conclusions

Five (Na, pH, OM, AN, SK) of the 10 soil indicators were selected as the MDS by PCA and external environmental factors, and the comprehensive quality index method was used to evaluate the soil quality in the Kenli District of the Yellow River delta. Based on the results of this study, the conclusions were as follows:

• (1). The average SQI of Kenli is 0.52, which is grade IV soil, and significant differences in the SQIs are noted among different land use types. The soil quality of agricultural land is greater than that of the natural land type as well as the bare land type. There is a certain law of the spatial distribution of soil SQI in the Kenli District, which mainly demonstrates that the closer the location is to the Yellow River, the better the soil quality will be. The soil quality of the eastern township is relatively greater than that of the western townships.

• (2). The analysis of the contribution of the MDS indicator to the SQI of each land type and township demonstrates that Na nitrogen and OM are the main indicators of limiting SQI in natural and bare land, and the main limiting factors of farmland soil quality are SK and pH. The limiting indicators of SQI vary among towns. Therefore, in the process of land development and utilization, each township should formulate and implement different management and development strategies according to its own characteristics and land types.

• (3). Within the study area, all NDVI, Ds, and Dr have different degrees of correlation with SQI. NDVI is positively correlated with SQI, and the two promote each other. Influenced by the freshwater and eutrophic fluvial marine sediments in the Yellow River, the closer the land is to the Yellow River and the coast, the better the soil quality is. In addition, no correlation was noted between elevation and soil quality in the study area.

Compared with similar studies in the Yellow River delta [29,31], this study describes in detail the spatial distribution of the soil quality of different grades in the county area, as well as the main limiting indicators of soil quality in towns and villages of different land types, and suggestions on the management and utilization of different land types in the township units of this region are proposed.

Author Contributions

Conceptualization, M.J., L.X. and X.C.; methodology, M.J. and H.Z.; software, H.F.; validation, L.X. and X.C.; formal analysis, M.J.; investigation, M.J. and X.C.; resources L.X. and X.C.; data curation, M.J. and H.Z.; writing—original draft preparation, M.J.; writing—review and editing, L.X.; visualization, H.F. and M.J.; supervision, L.X. and X.C.; project administration, L.X. and X.C.; funding acquisition, L.X. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research and Development Program of China (2018YFE0206400 and 2018YFC0407606), the National Natural Science Foundation of China (41971137 and 41771235), the Special Fund for Key R&D Program of Shandong (2017CXGC0316), the STS Key Projects of the Chinese Academy of Sciences (KFJ-STS-QYZD-098), the Science and Technology Planning Project of Qinghai Province (2019-HZ-818), the Nanjing Institute of Geography & Limnology (NIGLAS2019QD005) and the National Natural Science Foundation of China (U1806215).

Acknowledgments

We convey our sincere thanks to Kun Yan (Yantai Institute of Coastal Zone Research, CAS) and Yunpeng Sun (Institute of Soil Science, CAS). We also thank Congwen Wu for helping with the experiments.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Geographic locations of soil sampling sites with different vegetation types.

Figure 2. Flowchart of the study procedure.

Figure 3. Cluster–Pearson correlation coefficient heat map among soil properties. * Indicates significance level at p ≤ 0.05. ** Indicates significance level at p ≤ 0.01. N indicates no significance. The Pearson correlation analysis requires factors to conform to a normal distribution, so a natural logarithmic transformation is used for indicators that do not conform to a normal distribution.

Figure 4. The SQI value of different land use types (A) and different townships with the individual contributions of each MDS indicator (B). Error bars denote the standard deviation of the overall indicator value. In B, OM, AN, and SK are not significantly different, the Na contribution of cornfield is significantly greater than that of nudation, and the pH contribution of grassland is significantly greater than that of forest.

Figure 5. The SQI value of different townships with the individual contributions of each MDS indicator. Different lowercase letters represent values that are significantly different (p ≤ 0.05) between different productivity levels. Where (A) is the SQI value of different townships, and (B) is the individual contributions of each MDS indicator.

Figure 6. Omni-directional variograms calculated for SQI distribution.

Figure 7. Spatial distribution of soil quality in the Kenli District in 2013. The grades in the figure refer to the grade classifications of soil quality. According to the SQI > 0.7, 0.6–0.7, 0.5–0.6, 0.4 to 0.5 and ≤0.3, the soil quality can be divided into levels I to VI, and I is the best level in terms of soil quality. The interpolation graph is drawn by Golden Software Surfer 15 and color levels are assigned according to classification criteria.

Figure 8. Correlations between the SQI and NDVI (A), Dr (B), Ds (C), and Elev (D) of 127 sampling points in the Kenli District.

Table 1

The numbers of sampling points in different vegetation types and towns.

Vegetation Types	Huanghekou	Yong’an	Kenli	Xinglong	Shengtuo	Dongji	Haoji	Sum
Cash crop	2	0	1	0	2	0	0	5
Corn field	1	0	3	0	2	0	1	7
Cotton field	26	13	13	1	7	5	1	66
Forest	5	2	4	1	0	2	0	14
Grass	2	1	1	1	1	0	0	6
Nudation	2	1	1	0	2	0	0	6
Others	5	3	1	0	1	0	0	10
Paddy land	2	0	1	0	1	0	0	4
Reed	4	0	1	0	1	1	0	7
Suaeda salsa	2	0	0	0	0	0	0	2
Sum	51	20	26	3	17	8	2	127

Table 2

Descriptive statistics of the soil properties and external environmental factors.

Variables	Mean	SE	Min	Q1	Median	Q3	Max
NDVI	0.50	0.01	0.13	0.39	0.51	0.59	0.82
Elev	4.45	0.23	−0.94	2.71	4.50	6.17	12.37
Ds	23,276	1201.46	1040	12,204	20,743	30,374	55,683
Dr	7360	468.48	464	3080	6418	10,547	24,817
OM (g/kg)	10.15	0.38	1.19	6.85	10.17	12.48	21.77
pH	8.50	0.02	7.91	8.35	8.49	8.65	9.31
EC (mS/cm)	0.74	0.19	0.08	0.28	0.66	1.65	10.26
AN (mg/kg)	56.59	8.14	9.83	33.97	49.30	90.45	503.13
Cl (mg/kg)	462.76	359.57	18.76	79.76	520.97	1866.50	20,303.66
SO4 (mg/kg)	328.89	55.23	32.47	167.62	345.98	709.29	5634.49
Na (mg/kg)	339.35	153.57	20.51	104.30	359.18	1017.57	9378.79
Mg (mg/kg)	43.66	24.74	4.13	14.13	28.03	105.12	1746.39
Ca (mg/kg)	157.14	35.97	26.19	78.92	145.07	281.40	3237.80
SK (mg/kg)	15.99	3.68	2.55	9.53	15.48	22.43	258.96

(Where Mean is the average, SE is standard error, Min is minimum value, Max is maximum value, Q1 and Q3 are the upper and lower quartiles, respectively Because EC, AN, Cl⁻, Na⁺, SO₄²⁻, Ca²⁺, Mg²⁺, and SK do not obey a normal distribution, the mean is represented by geometric mean. Elev is the altitude of the sampling point, NDVI is normalized difference vegetation index, Dr and Ds represent the nearest straight-line distance from the sampling point to the Yellow River and the coastline, respectively, and the meaning of these abbreviations is the same in the tables and figures below).

Table 3

Results of PCA for soil quality indicators.

PCs	PC1	PC2	PC3
Eigenvalue	5.595	1.379	1.048
Percent	55.947	13.794	10.485
Cumulative percent	55.947	69.741	80.226
Eigenvectors
Cl	0.966	−0.144	0.036
SO4	0.817	0.013	−0.140
Na	0.909	−0.194	0.132
Mg	0.969	−0.053	−0.104
Ca	0.867	0.107	−0.261
EC 1:5	0.970	−0.103	0.028
pH	−0.495	−0.611	0.288
AN	0.055	0.661	0.560
SK	0.529	0.167	0.637
OM	−0.095	0.677	−0.358

Table 4

The coefficient of determination (R²) between the indicators and external environment.

Indicators	NDVI	Elev	Ds	Dr	OM	pH	lnEC	lnAN	lnCl	lnSO4	lnNa	lnMg	lnCa	lnK
NDVI	1	0.01	−0.08	−0.11	0.15	0.03	−0.21 *	0.14	−0.20 *	−0.09	−0.22 *	−0.15	−0.06	−0.21 *
Elev	0.01	1	0.75 **	−0.36 **	−0.14	0.17	−0.26 **	0.03	−0.30 **	−0.21 *	−0.31 **	−0.29 **	−0.24 **	−0.30 **
Ds	−0.08	0.75 **	1	−0.36 **	−0.14	0.24 **	−0.22 *	−0.11	−0.24 **	−0.16	−0.24 **	−0.32 **	−0.28 **	−0.42 **
Dr	−0.11	−0.36 **	−0.36 **	1	−0.05	−0.24 **	0.26 **	−0.08	0.27 **	0.16	0.26 **	0.31 **	0.27 **	0.25 **

(* Indicates significance level at p ≤ 0.05. ** Indicates significance level at p ≤ 0.01.).

Table 5

Indicators included in the MDS.

Group	Indicator	Norm	R2				Normal Transformation					Value	Included
			NDVI	Elev	Ds	Dr	Norm	NDVI	Elev	Ds	Dr
1	EC1:5	2.31	0.08	0.10	0.10	0.13	1.00	0.72	0.91	0.87	0.83	4.33	No
1	Cl	2.30	0.09	0.10	0.11	0.15	1.00	0.84	0.97	1.00	0.96	4.77	No
1	SO4	1.93	0	0.03	0	0	0.84	0	0.27	0	0	1.11	No
1	Na	2.18	0.11	0.11	0.11	0.15	0.94	1.00	1.00	0.97	1.00	4.91	Yes
1	Mg	2.30	0	0.08	0.10	0.11	1.00	0	0.71	0.90	0.69	3.29	No
1	Ca	2.08	0	0.04	0.05	0.05	0.90	0	0.39	0.48	0.35	2.12	No
2	OM	0.92	0	0	0	0	0.40	0	0	0	0	0.40	Yes
3	pH	1.43	0	0	0.06	0.06	0.62	0	0	0.52	0.38	1.14	Yes
4	AN	0.99	0	0	0	0	0.43	0	0	0	0	0.43	Yes
5	SK	1.43	0.08	0.03	0.07	0.05	0.62	0.69	0.31	0.62	0.33	2.24	Yes

Table 6

Standard scoring functions, communality and weight of soil quality indicators of the minimum data set (MDS).

Indicator	SSF	Value	Weight
Na	SSF3	0.881	0.241
pH	SSF3	0.690	0.189
OM	SSF1	0.623	0.170
AN	SSF1	0.772	0.211
SK	SSF1	0.691	0.189

Table 7

The soil quality comprehensive index (SQI) value and quantitative proportion of different land use types in each soil quality grade.

Land Use	Points	Means	CV	I (>0.7)		II (0.6–0.7)		III (0.5–0.6)		IV (0.4–0.5)		V (0.3–0.4)		VI (<0.3)
Total	127	0.52 ± 0.14	0.26	13	10.2%	25	19.7%	30	23.6%	37	29.1%	17	13.4%	5	3.9%
Cornfield	7	0.60 ± 0.2 a	0.33	3	42.9%	1	14.3%	0	0.0%	2	28.6%	0	0.0%	1	14.3%
Cotton	67	0.55 ± 0.14 a	0.26	9	13.4%	18	26.9%	17	25.4%	12	17.9%	9	13.4%	2	3.0%
Grassland	7	0.47 ± 0.08 ab	0.17	0	0.0%	1	14.3%	2	28.6%	3	42.9%	1	14.3%	0	0.0%
Forest	14	0.48 ± 0.07 ab	0.15	0	0.0%	1	7.1%	4	28.6%	8	57.1%	1	7.1%	0	0.0%
Reed	7	0.45 ± 0.14 ab	0.31	0	0.0%	1	14.3%	1	14.3%	2	28.6%	2	28.6%	1	14.3%
Nudation	6	0.36 ± 0.07 b	0.20	0	0.0%	0	0.0%	0	0.0%	2	33.3%	3	50.0%	1	16.7%

Notes: Means followed by the same letters in the rows are not different at p < 0.05 according to Tukey’s test.