Soil Quality and Evaluation of Spatial Variability in a Semi-Arid Ecosystem in a Region of the Southeastern Iberian Peninsula (Spain)

Abstract

In the last two decades, as the importance of soil has been recognized as a key component of any ecosystem, there has been an increased global demand to establish criteria for determining soil quality and to develop quantitative indices that can be used to classify and compare that quality in different places. The preliminary estimation of the attributes involved in soil quality was made taking into account the opinion of the experts and our own experience in a semi-arid ecosystem. In this study, 16 soil properties have been selected as potential indicators of soil quality, in a region between Campo de Montiel and Sierra de Alcaraz (Spain): sand and clay percentage, pH, electrical conductivity (EC), soil organic carbon (OC), extractables bases of change (Na, K, Ca and Mg), cationic exchange capacity (CEC), carbonate calcium equivalent (CCE), bulk density (BD), water retention at 33 kPa field capacity and 1500 kPa permanent wither point (GWC33 kPa and GWC1500 kPa), coefficient of linear extensibility (COLE) and factor of soil erodibility (K). The main objective has been to develop an adequate index to characterize the quality of the soils in a semi-arid Mediterranean ecosystem. The preliminary estimation of the attributes involved in soil quality was made considering the opinion of the experts and our own experience in semi-arid ecosystems. Two indicator selection approaches have been used to develop the Soil Quality Index (SQI) (total data set -TDS- and minimum data set -MDS-), scoring functions (linear -L- and nonlinear -NL-) and methods (additive -A-, additive weighted -W- and Nemoro -N-. The quality indices have been calculated, considering the properties of the soil control section (between 0 and 100 cm depth), using 185 samples, belonging to horizons A, B and C of 51 soil profiles. The results have shown that the election of the soil properties, both of the topsoil and subsoil, is an important help in establishing a good relationship between quality, soil functions and agricultural management. The Kriging method has been used to determinate the spatial distribution of the soil quality grades. The indices that best reflect the state of soil quality are the TDS-L-W and TDS-L-A should go as sub-indices, as they are the most accurate indices and provide the most consistent results. These indices are especially indicated when carrying out detailed or semi-detailed studies. However, the MDS-L-W and MDS-L-A should go as sub-indices, which use only a limited number of indicators, are best for large-scale studies. The indicators with the greatest influence on soil quality for different land uses and those developed on different rocks, using linear scoring functions, are the following: (Clay), (GWC1500 kPa) and (Ca). These results can also be expressed as follows: the best soils in this region are deep soils, with a clay texture, with high water retention and a neutral or slightly basic pH. However, the indicators with the greatest influence on soil quality, using nonlinear scoring functions, are: (OC Stock), (Ca) and (CaCO₃). In other words, the most important indicator is the organic carbon content, which is not logical in the case of a region in which the soils have an excessively low SOC content (0.86%).

1. Introduction

Soil has long been recognized as nonrenewable because, due to improper use or poor management, it can erode in a relatively short period of time, with little opportunity for regeneration. Even when the soil is contaminated, remediation may not be possible and returning it to its pre-contamination condition is impossible [1].

Over the past two decades, as the importance of the soil as a key component of natural and human-influenced ecosystems was recognized, the demand for establishing criteria for determining soil quality and developing quantitative indices that can be used to classify and compare the quality of soils in different places or in the same place over time [2]. It must be recognized that soils frequently perform several basic functions simultaneously: biomass production, the protection of human beings and the environment (the soil acts as a filter, buffer and transformation of harmful substances that can be spilled on it); the soil provides a habitat for numerous organisms and microorganisms (high biodiversity and genetic reserve); the soil constitutes the physical environment for the development of urban infrastructures; it is a source of raw materials and the soil acts as a carbon reserve. Soil functions are closely related to soil quality or, in other words, by indirectly evaluating soil quality we are evaluating soil functions. Any evaluation of the quality of soils must consider this multifunctional role. Soils provide a broad set of vital ecosystem services, but soils are under threat throughout the world. To avoid further degradation of soils and, consequently, the provision of soil-based ecosystem services it is essential to conserve ecosystem services, in this case soil-based services, because they are considered as the benefits that an ecosystem brings to society and that improve people’s health, economy and quality of life. When we are unable to conserve them, their degradation leads to significant damage to human well-being. Various soil-related ecosystem services can be considered: (a) provisioning, those referring to the amount of goods or raw materials that an ecosystem offers, such as wood, water, food or fibers; (b) regulatory, we can mention the control of soil erosion and climate change (soil as a carbon sink) and (c) support, such as the biodiversity of the soil ecosystem [3]. Any evaluation of soil quality should consider this multifunctional role. Within the European Union there is still no consensus on how to define soil quality. It is evident that at the present time there is a continuous demand for soil quality indices. Soil quality assessment has become an important activity in the last few decades, which is expected to increase in importance as we realize the need to protect and preserve the soil and its ability to maintain its functions. Given the complex nature of the soil and the exceptionally large number of soil properties that can be determined, it is important to be able to select the properties that are most appropriate for calculating these indices. It is logical that, depending on the nature of the basic function of the soil that is considered, the properties to select must vary [4,5].

Although many conceptual models have been proposed to assess soil quality, there is no universal method [6,7,8,9]. Linear and nonlinear scoring methods are the most widely used today [6,10]. Soil quality assessments have sometimes been performed using topsoil properties [10,11,12] and studies using complete soil profile data are limited [13,14]. The properties of the topsoil can be easy to measure and evaluate, but they provide incomplete information, as soil functions are driven by pedogenic processes in the “control section” of the soil [15].

In a previous study [16] different soil quality indices have been calculated, in a comparative way, not only using the properties of the arable layer (topsoil, between 0–25 cm), but also the properties of all the horizons of the soil profile (between 0 and 100 cm). Evaluation of the soil quality index using only the properties of the soil surface horizon (0–25 cm) provides incomplete information, since the functions of the soil and its productivity are influenced by the properties of both the surface and the subsurface horizons of the soil. Logically, the presence of a lithic contact or a petrocalcic horizon below horizon A is not the same as the presence of a horizon B with clay texture and great thickness.

For example, in a Mediterranean semi-arid ecosystem with a predominance of a soil association made up of Entisols and Alfisols. The Entisols, thin soils with an A horizon developed on compact and hard rock, the presence of this lithic contact prevents the root development of the crops and that soil must have low fertility. On the other hand, the Alfisoles, deep soils, with a balanced texture in the superficial horizon, have a subsurface horizon with a clayey texture, high cation exchange capacity and high water retention (in a region where soils have a xeric moisture regime, which implies a deficit in the water balance in the soil between the months from April to October). Logically, on this type of deep soil the crops will not have root development problems and the storage of water in the soil, in the root zone, during the wet season is essential to maintain the restoration of vegetation and growth of the crop during the dry season [17]. For this reason, in this work we have studied the soil quality index using the properties of all the horizons of the soil profile (between 0 and 100 cm).

Numerous soil quality assessment studies focus on physical and chemical indicators of the soil and are rarely described using biological indicators [18,19,20]. In the present study, there is no database on the biological activity of the soils in the studied region.

Frequent recent studies have reported a correlation between the biological and physicochemical properties of the soil (organic matter, pH, nutrient content, removable capacity for change, moisture retention, etc.) [21,22]. That is, some biological properties of the soil can be explained from the aforementioned physicochemical properties, which have a very important influence not only on the biological activity of the soil, but also on the growth of plants, on the nutrient cycle available for plant absorption and in improving soil fertility [23].

The main objectives of this study were: (1) to determine the appropriate indicators to assess the soil quality of this region, focusing on the support function of natural ecosystems and food and fiber production and using a statistical approach and linear relationships; (2) to evaluate the quality of the soils in a semi-arid ecosystem in Spain, using two indicator selection methods (TDS and MDS), two types of equations (linear and nonlinear) and three methods (Weighted additive, Additive and Nemoro). The quality indices have been calculated, taking into account the properties of the soil control section (between 0 and 100 cm depth), using 185 samples, belonging to horizons A, B and C of 51 soil profiles; (3) to establish a relationship between the quality index obtained and the land uses (agricultural, grassland and forestry) and geological material from which the soil is formed (quartzites–slates, sandstones, limestones–dolomites and clays–marl); (4) spatial predictions of soil quality index in order to identify the worst-quality areas and avoid further degradation with their agricultural use.

2. Materials and Methods

2.1. Description of the Study Area

The study area (Figure 1) is located in the surroundings of the limit that separates the provinces of Ciudad Real, Albacete and Jaén (Spain). The coordinates of the studied area are as follows: East–West-4,261,000–4,292,000; North–South-513,000–550,000; UTM Zone-30 y Datum-ETRS89. Geographically, the region studied is located at the confluence of Campo de Montiel, Sierra del Relumbrar and Sierra de Alcaraz. It covers a total area of 810 km². The highest level of the studied sector (1795 m.a.s.l.) is located in the Sierra de Alcaraz. The maximum level of the Sierra del Relumbrar is the Pilas Verdes peak (1151 m.s.n.m.). The minimum altitude is at the confluence of the Villanueva de la Fuente and Guadalmena rivers, at a level of 661 m.s.n.m.

The region object of the present study is located in three different geological units (Figure 2):

• The Central Iberian Zone, belonging to the Hercinian chain of the Iberian Massif (Sierra del Relumbrar): according to the division adopted in the Tectonic Map of the Iberian Peninsula and the Balearic Islands [24]. The Lower Ordovician is represented by a quartzite, with a thickness of a few hundred meters. The Middle Ordovician is made up, everywhere, of dark slates, a few hundred meters thick, with fauna from the Llanvirn and the Llandeilo. The highest parts of the Llandeilo and the Caradoc are formed by sandstones and shales. Regarding the structure, the most striking element is formed by the long, shapeless nuclei that run longitudinally in the area. They are nonconformal with an axial surface quite close to the vertical, although towards the S and SE the vergences tend towards the S.
• Mesozoic (and paleogen) covert of the Iberian Massif, little or no deformation (Campo de Montiel): the paleozoic materials of the Sierra del Relumbrar constitute the plinth on which mesozoic and tertiary materials have been deposited. The plinth, in numerous points, is hidden under these discordant materials that constitute an unfolded platform cover. It is a flat area, with a tabular relief, very extensive and with a dolomitic substrate. On both sides of the Sierra del Relumbrar are two broad valleys modeled on soft materials from the Triassic age (clays, sandstones, marls and plasters). The aforementioned Sierra del Relumbrar (with altitudes of approximately 1100 m) emerges in the center of two broad valleys (with altitudes of 800 to 850 m).
• The Prebética Zone (Sierra de Alcaraz): the Cordillera Bética constitutes the great structural unit of the South of the Iberian Peninsula. It is part of the set of alpine ridges of the Western Mediterranean. Three main bands or units are distinguished, lengthened in the WSW to ENE direction: (1) the Prebética Zone; (2) the Subbética Zone and (3) the Bética Zone.

The pre-Mesozoic plinth does not emerge in the Prebética Zone. It includes the Mesozoic and Tertiary terrains [25,26] up to the Lower Miocene inclusive. The facies are continental and shallow marine waters. The only materials that emerge within the sector studied in this area belong to the Lías and possibly the Dogger. The Lower Lias is made up of limestones and dolomites and, sometimes, levels of evaporites. In the Middle–Upper Lees there are levels of evaporites and abundant clays. The Dogger is made up of limestone and dolomites. Scale tectonics is one of the main characteristics of Prebético [27].

The climate of this region is characterized by long and cold winters and short, hot and dry summers, with an average rainfall of 663 mm per year and an average annual temperature of 14 °C. Precipitation is mainly concentrated in the winter period (December). These soils have a mesic temperature regime and a xeric moisture regime.

On the soils of the study area, classified according to the Soil Taxonomy [28] and World Reference Base for Soil Resources of FAO [29], Alfisols (Luvisols), Inceptisols (Cambisols), Entisols (Fluvisols, Regosols and Leptosols), Mollisols (Phaeozems, Leptosols rendzic and Gleysols mollic) and Vertisols (Vertisoles) predominate.

To carry out this study, a series of auxiliary data (Figure 3) was used, which correspond to factors of considerable influence within the selection process of the indicators that make up the soil quality index, among which are the map land use (the area is characterized by the coexistence of areas with natural vegetation, little or nothing degraded, with others were dedicated to intensive agriculture, and were moderately degraded). The slope map and terrain elevation have been made from the 5 m resolution Digital Terrain Model (DEM) of the studied area, using the ArcGIS 10.5 software package. Slopes were classified into five classes, that is, <5° (light slope), 5–12° (light slope), 12–20° (moderate slope), 20–29° (steep slope) and very steep slope (>29°). The terrain elevation was also classified into several classes, with an altitude between 661 and 1795 m. The soil erodibility map (K factor) shows the lowest erosion rates (K = 0.23–0.41) in soils with high useful depth (Alfisols), located on flat slopes and developed on quartzite debris or dolomites, with a high content of stony and/or rocky fragments, which act as a protective shield against erosion. The highest rates of erodibility (K = 0.51–0.81) occur in soils of small to moderate development (Entisols–Inceptisols), with little useful depth, located on steep slopes and developed on dolomites, slates, quartzites or sandstones, with little or no stoniness.

On the soils of the studied area a continuous vegetal carpet develops that occupies most of the area. Approximately 50% of the vegetation is autochthonous (cleared holm oak forest, scrubland and grasslands “Mediterranean climax forest”), while the other half is made up of man-made crops (cereals and olive groves). The need to feed more and more abundant settlers, logging—sometimes indiscriminate—and the subsequent erosion have reduced the primitive vegetation to what can be seen today on the land use map (Figure 3). Small remains of the original holm oak are still preserved, and where it has been modified, the successive stages of degradation appear, ranging from a thicket to a pasture. One of the stages of degradation of the primitive oak forest gives rise to two different types of scrub: the jarales (sticky rockrose, steppe, cantueso, etc.) that settle on slates and quartzites, while the kermes oak (kermes oak, rosemary, gorse, lavender, etc.) grows on limestone, loamy and dolomitic soils. The grasslands are formed by low-growing plant formations that are arranged forming more or less dispersed bushes depending on the nature of the terrain and its degree of erosion (thyme, helianthmus, grasses, etc.). In the Sierra de Alcaraz, on the SE border of the studied area, a small native pine forest (larice and resin pine) grows on limestone and dolomites.

2.2. Soil Sampling and Analysis

For the general characterization of the soils of the studied area, 185 samples were taken, belonging to horizons A, B and C from 51 soil profiles, that is, an average of three to four samples have been collected in each profile.

The samples were air-dried, minced and sieved through a 2 mm sieve, before performing chemical and physical analyses (

Table 1. Methods used in laboratory analyses for the selected indicators.

Indicator	Method
Granulometric analysis (% of sand and clay)	Robinson’s pipette [30]
Soil organic carbon (OC)	Dichromate oxidation [31]
Bulk density (BD)	Paraffin [32]
Water retention at 33 and 1500 kPa (GWC)	Pressure membrane [33]
Cation Exchange capacity (CEC)	Ammonium acetate at pH 7 [30]
Removable exchange bases (Na, K, Ca and Mg)	Ammonium acetate at pH 7 [30]
Electrical conductivity (EC)	Saturated soil paste and conductivity meter [34]
pH	Potentiometric method (1:1-soil-water-)
CaCO3 equivalent (CCE)	Bernard calcimeter [30]
Coefficient of Linear Extensibility (COLE)	Richards membrane [30]
Soil erodibility (K)	Equation of Wischmeier and Smith [35]

). The soil erodibility factor (K) has been calculated using the Equation (1):

100 K = [10 − 4 × 2.71 × T1.14 × (12 − OM)] + 4.2 × (E − 2) + 3.23 × (P − 3)]

(1)

where T is related to the texture factor in the shallowest 15 cm: T = [(100 − Ac) × (L + Armf)], with Ac = Clay, L = Slime, Armf = very fine sand, OM = organic matter (OC × 1.72), E = structure parameter and P = permeability.

Soil organic carbon (OC) is a key indicator of soil quality [36] and of the general productivity of the soil [37], and increases cation exchange capacity, aggregation and water retention [38].

A basic function of the soil is to act as a carbon sink. OC stored in soils represents the largest reserve of terrestrial carbon. Arable soils generally have low soil organic carbon values, while values are highest under permanent groundcover. Conversion of natural to agricultural land is estimated to result in 50–100 Pg losses of soil organic carbon worldwide over the past 200 years [39].

The soil organic carbon stock has been calculated as the product of three variables, SOC, BD and stoniness, using the following equation (Equation (2)):

Stock of OC = OC × BD × D × (1 − S)

(2)

where OC is the percentage of organic carbon (g 100⁻¹ g), BD is the bulk density (g cm⁻³), D is the thickness of the surface layer (m) and S is the percentage of the volume of the horizon occupied by thick shards (g 100⁻¹ g).

The sixteen selected indicators, which have been included in a TDS, have been chosen because they have a strong influence on soil quality and have been suggested by numerous authors due to their influence on soil fertility, nutrient supply, root growth, soil porosity, etc. [9,40,41,42]. In addition, the soil erodibility factor (K) has been included as part of the quality index in order to also consider human impacts due to agricultural practices and land use change [43]. In addition, the meteorological phenomenon known as “cold drop”, which is associated with extremely violent downpours and storms, typical of Spain, is generally associated with erosion processes that lead to significant soil loss, causing numerous adverse and long-lasting effects on soil properties [44].

Initially, the use of many different soil properties to obtain a quality index can greatly complicate the interpretation and synthesis of the results. However, experience indicates that the combination of these properties in a single general index makes the evaluation more meaningful and practical [10,11].

2.3. Statistical Analysis

For the determination of the different routine parameters of statistical interest, SPSS v.25 software was used. In addition, factor analysis (FA) was carried and to compare the different indices, the precision of the classification for each quality grade (very high, high, moderate, low and very low quality) was evaluated using the Kappa statistic and the correlation coefficients [45,46].

Geostatistical analysis allowed mapping of the spatial distribution of the soil quality classes identified using spatial interpolation methods.

The GIS analysis-ArcGis v10.9 with extensions Spatial analysis tools and geostatistical analyst, used to calculate the degree of spatial variability of the different calculated quality indices, has been calculated using “kriging”.

2.4. Score of the Indicators

During the scoring of the indicators, data normalization is required—scores ranging from 0 to 1 [6,10,47], since the indicators are generally expressed with different numerical scales. Liu et al. [48] used the linear scoring method (L) to normalize the data. This method establishes the linear relationship between the quality score and the measured data based on the indicator’s sensitivity to changes in soil quality. Based on the aforementioned sensitivity of the indicator in soil quality, three types of functions were applied:

4.

To CEC, GWC33 kPa, GWC1500 kPa, content of CaCO₃, extractable bases (Na, K, Ca and Mg) and OC Stock parameters, due to their role in soil fertility and water availability [6,42,49], the “More is better” function was applied. In this case, each indicator value was divided by the highest value, so that the highest value received a score of 1; that is, the following linear scoring curve was used (Equation (3)):

SL = X/Xmax

(3)

where SL is the linear score of the soil indicator, X is the measured value of the soil indicator and Xmax is the maximum value of the soil indicator.

5.

To the K factor “erosion”, EC and BD parameters, related with soil porosity and soil degradation [10,50], the “Less is better” function was used. The lowest value was divided by each indicator value so that the lowest value received a score of 1; that is, the following linear scoring curve was used (Equation (4)):

SL = Xmin/X

(4)

where SL is the linear score of the soil indicator, X is the measured value of the soil indicator and Xmin is the minimum value of the soil indicator.

6.

To the clay and sand content, COLE and pH, the “Optimal range” function was applied. In this case, the threshold values or the optimal ranges were identified: 40 and 80%, respectively, for the clay and sand contents; 0.060 for COLE and 6.5–7.5 for pH [45,47,51]. Scores were assigned using the more is better or the less better function, depending on whether the indicator value was below or above the optimal range.

Other studies did not find a linear relationship between quality scores and indicator values, and therefore developed the nonlinear scoring (NL) method to normalize the data [51,52]. For nonlinear scoring, the following sigmoidal curve (Equation (5)) was used as follows:

SNL = 1/1 + (X/Xm)²

(5)

where SNL is the nonlinear score of the soil indicator, X is the measured value of the soil indicator, Xm is the mean value of the soil indicator and b is the slope of the equation and is set to −2.5 for a “more is better” and 2.5 for a curve “less is better” [6,52].

To determinate the soil quality index, the weighted average of all the horizons of each soil profile studied (between 0 and 100 cm depth) was calculated. In this way, a single value was obtained for each index.

2.5. Selection of Minimun Data Set

The indicator selection methods, Total Data Set (TDS) and Minimum Data Set (MDS) are mostly used in the evaluation of soil quality.

The results of this study have identified a better estimate of the quality of the soil by applying the SQI-TDS index, in which all the properties of the soil are used and, therefore, it takes a long time and high economic costs to carry out the analysis from the laboratory. This index is especially indicated when carrying out detailed or semi-detailed studies. However, the SQI-MDS index, which uses only a limited number of indicators, reduces the cost of the analysis and, in turn, increases the sampling density.

The MDS method is recommended for evaluating soil quality in large-scale studies and in developing countries, where the measurement of indicators should be carried out as economically as possible and with minimal infrastructure.

Selecting a MDS is detrimental to the loss of information from the indicators that are not selected, but avoids problems such as information redundancy and tedious laboratory work [53].

Andrews et al. [10] and Imaz et al. [54] obtained a minimum set of indicator data from a total set of data using factor analysis, which provided high consistency in the evaluation of soil quality. During the FA, the Varimax rotation method was used in order to obtain a solution, as simple as possible, of the matrix of “loadings” with which each variable contributes to each of the factors obtained. With this rotation, the variance of the loadings within each factor is maximized and, in addition, the loadings will tend to take high or low values and, simultaneously, each variable will tend to have high loadings in only one factor.

The number of factors has been selected such that the eigenvalues are greater than one or very close to one and the explained variance is greater than 71% [11]. During this process it has to be accepted that soil variables with high factor loadings are the soil properties that best represent changes in soil quality. In particular, soil properties with absolute values around 20% of the highest factor loading were chosen [10,55]. Therefore, a model with four factors was chosen.

2.6. Soil Quality Indices

After the qualification of the indicators, the scores of the selected indicators are combined into a soil quality index. The method used to calculate the soil quality index (SQI), considers the importance of each indicator and specifies the weight value of each indicator [50,56], with weight value assignments based on expert opinion or statistical analysis [46].

For each of the TDS and MDS methods, weighting values were assigned considering the communality of each indicator, obtained from factor analysis (FA). The communality value is representative of the portion of the variance explained by each indicator. This value varies between zero and one, and a high communality value leads to a higher contribution of the indicator to soil quality [57]. The weight values for each indicator considered in this study were calculated from the ratio of communality of each indicator to the sum of all indicators of communalities [46,58].

The indicator scores were integrated into the indices using additive methods (Equation (6)), weighted additives (Equation (7)) [46,48] and the Nemoro quality calculation method (Equation (8)), which evaluates soil quality based on the minimum and average scores of the indicators [45,59]:

$${\mathrm{SQI}}_{\mathrm{A}}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{Si}/\mathrm{n}$$

(6)

$${\mathrm{SQI}}_{\mathrm{W}}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{Wi}\mathrm{Si}$$

(7)

$${\mathrm{SQI}}_{\mathrm{N}}=\frac{\sqrt{{\mathrm{Pave}}^2+{\mathrm{Pmin}}^2}}{2}\times \frac{\mathrm{n}-1}{\mathrm{n}}$$

(8)

In these equations, SQI is the soil quality index, S_i is the indicator score (linear or nonlinear), “n” is the number of soil indicators in the total data set or the minimum data set, W_i is the weight of each indicator [10,46], P_ave is the average of the indicators selected at each sample point and P_min is the minimum of the scores of the indicators selected at each site.

The quality indices (SQI_A, SQI_W and SQI_N) were calculated in all the “n” indicators qualified and weighted in the TDS and MDS methods for each sample.

Combining indicator selection methods and scoring function approaches, twelve indices were compared in the present study (Figure 4).

For each soil quality index, five classes or grades (very low, low, moderate, high and very high) were calculated as follows: the rank of each quality index was divided by the desired number of intervals (five), and the result was used as the width of each interval. Adding this value to the lowest value of the corresponding index gave the upper limit of the first interval, and so on, until the highest range of the quality index was reached.

3. Results and Discussions

3.1. Physical and Chemical Properties

The selection of indicators to calculate the soil quality indices has been carried out with all the physicochemical properties that are generally measured in a regional soil study.

Table 2. Descriptive statistics of soil properties.

Properties	Average	Minimum	Maximum	Standard Deviation	Kurtosis	Coef. of Variation (%)
Sand %	31.91	2.00	96.00	18.72	0.42	58.68
Clay %	34.65	2.20	85.00	16.20	−0.34	46.76
OC %	0.86	0.02	3.75	0.82	2.31	96.06
Stock OC (mg C ha−1)	40.60	10.75	94.50	20.39	0.12	50.22
CaCO3 (%)	18.39	0.00	88.00	23.24	0.07	126.41
BD (g cm−3)	1.49	0.83	1.96	0.16	2.31	10.62
GWC33kPa (%)	27.17	6.00	43.50	7.77	−0.44	28.59
GWC1500kPa (%)	15.93	1.90	31.50	6.80	−0.63	42.66
COLE (cm/cm)	0.032	0.000	0.116	0.027	0.27	83.13
pH	7.32	4.50	8.60	0.88	0.48	11.99
Na (cmol kg−1)	0.44	0.00	5.60	0.59	35.35	133.36
K (cmol kg−1)	0.63	0.00	2.50	0.48	1.60	77.46
Ca (cmol kg−1)	28.26	0.30	65.30	16.69	−0.97	59.06
Mg (cmol kg−1)	3.21	0.00	15.20	2.96	2.09	92.10
CEC (cmol kg−1)	18.79	0.90	65.20	10.79	0.99	57.46
CE (ds m−1)	0.69	0.13	4.02	0.54	16.42	77.61
Factor K (t m2 h/ha J cm)	0.46	0.23	0.81	0.11	0.61	25.05

shows the main statistical values of the 16 indicators measured at each sampling point.

Most soil quality assessment studies focus on soil physical and chemical indicators and are rarely described by biological indicators [21,23]. Biological indicators such as activity dehydrogenase, β-glucosidase, urease, acid phosphatase, microbial biomass C, microbial biomass N and breathing rate are rarely used in regional soil studies, because the measurement of these properties in numerous samples is usually very expensive. In the present study, no biological property has been used because there is no database on the biological activity of soils in the studied region.

Frequent recent studies have reported a correlation between the biological and physicochemical properties of the soil (organic matter, pH, nutrient content, extractable exchange capacity, moisture retention, etc.) [24,25]. That is, some biological properties of the soil can be explained from the aforementioned physical–chemical properties, which have a very important influence not only on the biological activity of the soil, but on the growth of plants in the nutrient cycle available for plant uptake, and in improving soil fertility [26].

3.2. Factorial Analysis

The results obtained from the factor analysis (FA) provided a solution with four factors, when considering the criterion of eigenvalues > 1. Once these factors were rotated, the loads were obtained (

Table 3. Reordered matrix of rotated factors.

	F1	F2	F3	F4
Eigenvalues	6.047	2.484	1.642	1.189
Variance (%)	37.793	15.524	10.260	7.430
Cumulative variance (%)	37.793	53.317	63.577	71.007
Soil properties
GWC1500kPa (%)	0.928
Clay (%)	0.916
GWC33kPa (%)	0.912
Sand (%)	−0.890
CEC (cmol kg−1)	0.870
COLE (cm/cm)	0.821
CaCO3 (%)		0.837
Ca (cmol kg−1)		0.823
pH		0.799
Na (cmol kg−1)		0.592
Stock OC (mg C ha−1)			0.712
Factor “K” (t m2 h/ha J cm)			−0.689
BD (g cm−3)			−0.525
Mg (cmol kg−1)				0.843
CE (ds m−1)				0.819

In bold: factor loads selected as MDS.

). The Kaiser–Meyer–Olkin (KMO) test resulted in a value of 0.710, which indicates an excellent adequacy of the sampling carried out. Then, the description of the values found was made.

Factor 1 has high loads in % clay and % sand and in those that would cavate with them, such as water retention parameters, exchange capacity and COLE. The sand is presented with negative charges. This factor is clearly granulometric. This factor is the one with the greatest weight in the identity of the soils, since it explains 37.79% of the variance.

Factor 2 is “compositional” and involves the variables that express the content of calcium carbonate, Ca and extractable Na. A variable that covariates with the carbonates also appears, such as the pH (which increases with the increase in the content of calcium carbonate, Ca and extractable Na).

Factor 3 only carries a high load on the variable organic carbon stock. It should be noted that the apparent density has a charge of the order of 0.525 (negative). This factor could be defined as “the increase in organic carbon leads to a decrease in the apparent density in the soils”. The soil erodibility factor “K” has a 0.689 (negative) charge.

Factor 4 carries high loads on the variables of extractable Mg and on the conductivity of the saturation extract. It can be affirmed that this factor expresses the salinity of the soils.

Considering the community analysis and the weight value of the indicators (

Table 4. Result of the estimated community and the value of the weight of each indicator in both TDS and MDS methods.

	TDS		MDS
	Communality	Weight	Communality	Weight
Sand%	0.803	0.071
Clay%	0.848	0.075	0.860	0.212
Stock OC mg C ha−1	0.570	0.050	0.108	0.027
Ca CO3%	0.732	0.064	0.872	0.215
BD g cm−3	0.354	0.031
GWC33kPa%	0.907	0.080
GWC1500kPa%	0.917	0.081	0.922	0.228
COLE cm cm−1	0.774	0.068
pH	0.708	0.062
Na cmol kg−1	0.581	0.051
K cmol kg−1	0.321	0.028
Ca cmol kg−1	0.899	0.079	0.891	0.220
Mg cmol kg−1	0.831	0.073	0.395	0.098
C.E.C. cmol kg−1	0.867	0.076
CE ds m−1	0.710	0.062
Factor “K” t m2 h/ha J cm	0.539	0.047

), the water retention at 1500 and 33 kPa, as well as the extractable Ca, the CEC and the clay content, received the highest weights (between 0.081 and 0.075); on the contrary, the extractable K, the apparent density and the erodibility factor “K” showed the lowest weights (between 0.028 and 0.047).

The minimum set of variables (MDS) that best describes the soils of the region should be a selection made up of one, two or three representatives of each factor with the highest loads. For this reason, the following indicators were selected: the selected soil parameters of F1 were GWC_1500kPa% and Clay%; the selected parameters of F2 were CaCO₃ and extractable Ca; the organic carbon stock was selected from F3 and the extractable Mg from F4.

3.3. Soil Quality Grades and Their Spatial Distribution

Each of the twelve soil quality indices obtained from the properties of all soil horizons (0–100 cm) classified the soil quality of the study area into five classes (Table 4 and

Table 5. Classification of grades of soil quality.

Index	Soil Quality Grade
	I (Very Low)	II (Low)	III (Moderate)	IV (High)	V (Very High)
SQITDS-L-W	0.31–0.37	0.37–0.43	0.43–0.50	0.50–0.56	0.56–0.63
SQITDS-L-A	0.33–0.38	0.38–0.44	0.44–0.49	0.49–0.55	0.55–0.61
SQITDS-L-N	0.22–0.26	0.26–0.30	0.30–0.33	0.33–0.37	0.37–0.42
SQITDS-NL-W	0.14–0.24	0.24–0.33	0.33–0.43	0.43–0.52	0.52–0.63
SQITDS-NL-A	0.17–0.26	0.26–0.35	0.35–0.43	0.43–0.52	0.52–0.62
SQITDS-NL-N	0.11–0.18	0.18–0.25	0.25–0.33	0.33–0.40	0.40–0.48
SQIMDS-L-W	0.08–0.21	0.21–0.33	0.33–0.46	0.46–0.58	0.58–0.72
SQIMDS-L-A	0.08–0.20	0.20–0.32	0.32–0.43	0.43–0.55	0.55–0.68
SQIMDS-L-N	0.05–0.13	0.13–0.21	0.21–0.29	0.29–0.37	0.37–0.46
SQIMDS-NL-W	0.04–0.18	0.18–0.31	0.31–0.45	0.45–0.58	0.58–0.73
SQIMDS-NL-A	0.11–0.23	0.23–0.35	0.35–0.46	0.46–0.58	0.58–0.71
SQIMDS-NL-N	0.07–0.16	0.16–0.25	0.25–0.33	0.33–0.42	0.42–0.52

).

Current digital soil mapping (DSM) techniques take advantage of advances in computer hardware, geographic information systems and statistical techniques. Geostatistical analysis, using ArcGIS 10.5 software, has made it possible to make a map, showing the spatial distribution of the different degrees of soil quality, using the kriging method of spatial interpolation. Interpolation is the process of predicting values to unknown sites, considering the information on the geographical location of the points, actually sampled [60].

A DSM approach has many advantages over conventional soil mapping approaches; For example, by leveraging increasingly free and easily accessible geospatial data sets, and in conjunction with predictive modeling techniques, soil map development can be automated to develop products that are more accurate than conventional soil maps, which can rarely be updated, due to the excessive time and cost involved. In addition, digital soil mapping (DSM) techniques include a set of useful tools that facilitate large-scale soil mapping in data-poor regions.

In Figure 5 it can be seen with the naked eye that the spatial patterns of soil quality derived from the 12 methods used are similar. The parent material played an important role in the spatial distribution of soil quality. A similar pattern can easily be seen when looking at both maps: the soil quality map and the geological one.

The quality of the soil in the study area is preferably moderate (Grade III) -green areas-, when using the TDS-L and MDS-L models. The maps made using the TDS-NL and MDS-NL models show a predominance of both “green areas” or “moderate quality” (Grade III) and “light blue areas” or “high quality” (Grade IV). Only a small percentage of the surface of the studied area, in all the models, has a very high-quality grade (Grade IV) “dark blue areas” and the areas that have soils of very low quality (Grade I) “red areas” (although in some maps the surface occupied by low-quality soils is 0.00 hectares, there are actually one to five soil profiles that are considered within this grade).

In the spatial distribution map, grade III (Moderate) of the TDS-L-W and TDS-L-A indices are the ones with the highest representation (442 and 498 km², respectively). However, the same grade III, in the TDS-L-N index only occupies 287 km² (

Table 6. Area, in km², corresponding to each grade of soil quality.

Data Set	Equation	Method	Soil Quality Grade
			I (Very Low)	II (Low)	III (Moderate)	IV (High)	V (Very High)
TDS	Linear	Weighted additive	0.19	151.72	442.38	186.36	31.02
		Additive	10.20	127.75	497.83	151.81	24.24
		Nemoro	15.36	260.09	286.70	221.87	27.57
	Nonlinear	Weighted additive	0.00	10.99	279.37	469.58	51.66
		Additive	0.00	12.50	239.42	497.28	62.40
		Nemoro	0.00	39.73	503.23	252.50	16.13
MDS	Linear	Weighted additive	0.00	84.68	392.10	294.66	40.15
		Additive	0.00	80.85	397.70	298.17	34.88
		Nemoro	0.00	125.27	492.48	172.71	20.97
	Nonlinear	Weighted additive	0.00	67.86	413.21	310.11	20.42
		Additive	0.00	59.11	383.52	360.02	8.94
		Nemoro	0.01	158.98	553.05	99.35	0.20

). The studied region occupies a total area of 810 km², therefore, grade III, in each of the TDS-L-W and TDS-L-A indices, occupies 54.5 and 61.5%, respectively. However, the same grade III, in the TDS-L-N index, only occupies 35.4% of the total extension.

In the indices made with a nonlinear equation (TDS-NL-W and TDS-NL-A), grade IV (High) is the one with the greatest extension on the map (470 and 497 km², respectively). However, the same grade IV, in the TDS-LN-N index only occupies 252 km². Grade IV, in each of the TDS-L-W and TDS-L-A indices, occupies 58.0 and 61.3%, respectively. However, the same grade IV, in the TDS-L-N index only occupies 31.1% of the total extension.

In the MDS-LW and MDS-LA indices, grade III (Moderate) is again the one with the largest area on the map (392 and 398 km²; that is, 48.4 and 49.1% of the surface total, respectively). However, the same grade III, in the MDS-L-N index occupies a very high area (492 km²), 60.7% of the total area.

In the indices made with a nonlinear equation (MDS-NL-W and MDS-NL-A), grade III (Moderate) is the one with the greatest extension on the map (413 and 383 km²; that is, 51, 0 and 47.3% of the total area, respectively). However, the same grade III, in the MDS-LN-N index occupies a very high area (553 km²), 68.3% of the total area.

Summarizing, the surface percentages of the different degrees of soil quality obtained in the TDS and MDS models are similar with the use of the Weighted Additive (W) and Additive (A) methods; however, with the Nemoro (N) method the results are inferior (they are underestimated).

In general, the studied area is an important area for the agricultural production of rainfed cereals (preferably wheat and barley), due to its high soil fertility, aptitude for agriculture and potential productivity. If we make a geographical distribution of the quality of the soil types in the studied region, three zones can be distinguished: (1) zone of very low–low quality (red and yellow colors), which corresponds to Entisols developed on quartzites, slates and sandstones, located in the Sierra del Relumbrar, in the valleys on Triassic sandstones that surround the Sierra del Relumbrar and on dolomites in the Mesozoic Covert of the Iberian Massif; (2) of moderate quality (green zone), located in the central part of the map, which closely corresponds to the predominance of the Inceptisols and Alfisols developed over dolomites and maciños of the Sierra de Alcaraz and over deposits of the glacis surrounding Sierra del Relumbrar and (3) high–very high quality (light blue and dark blue colors), which closely corresponds to the predominance of Inceptisols, Alfisols, Mollisols and Vertisols’ developed deposits of glacis and clays and marls from the Mesozoic Covert. and the Sierra de Alcaraz.

Considering the expert opinion of agricultural technicians and farmers, the soils of this region are considered predominantly of moderate–high quality, especially for rain-fed cereal cultivation. Digital soil maps, which identify areas with different soil quality grades, can be used to facilitate better land management and avoid soil degradation processes.

3.4. Validation of Quality Indices

To compare the different indices, the precision of the classification for each quality grade was evaluated using the Kappa statistic and the correlation coefficients [52,57]. For Kappa analysis, a value was calculated to show the following levels of agreement [61,62]: (1) null: <0.05; (2) very low: 0.05–0.2; (3) low: 0.2–0.4; (4) moderate: 0.4–0.55; (5) good: 0.55–0.7; (6) very good: 0.7–0.85; (7) near perfect: 0.85–0.99; and (8) perfect: 1.

The evaluation of the agreement of the degrees of quality determined by the TDS and MDS indices and between linear scoring methods: that is, between the TDS-L-W and MDS-L-W and between TDS-L-A and MDS-L-A, in both cases the same Kappa value (0.45 (moderate)). The evaluation of the agreement between the quality grades determined by the TDS-L-N and MDS-L-N indices has resulted in a Kappa value of (0.42 (moderate)).

When comparing the agreement between the nonlinear scoring methods, the Kappa statistical values are as follows: between the TDS-NL-W and MDS-NL-W indices, a Kappa value of (0.38 (low)) has resulted. Between the TDS-NL-A and MDS-NL-A indices a Kappa value of (0.40 (between low and moderate)) has resulted. Finally, between the TDS-NL-N and MDS-NL-N indices, a Kappa value of (0.29 (low)) has resulted.

According to the Kappa analysis, the lowest levels of agreement are presented in the indices calculated using a nonlinear score and the Nemoro method.

When performing the validation of the quality indices, using the TDS and MDS approaches (Figure 6) through the use of correlation coefficients, they have been highly correlated with each other (all R values are very high, between 0.81 and 0.85). The highest coefficient (R = 0.85) resulted when using the TDS and MDS-NL-A (additive and nonlinear scoring) models. When using the TDS and MDS-L-A models, and TDS and MDS-L-W (Additive and Additive weighted and linear scoring) are used, in both cases a coefficient very similar to the previous one (R = 0.83) results.

The lowest coefficient (R = 0.81) results from the use of the TDS and MDS-L-N models (Nemoro and linear score).

R values between 0.54 and 0.95 have been obtained by correlating the quality indices using linear and nonlinear scoring methods (Figure 7). Of the six correlations made, five have values between 0.54 and 0.68.

The highest coefficient (R = 0.95) has resulted when using the MDS-L-W and MDS-NL-W models (weighted additive models and linear and nonlinear scores, respectively).

The results obtained, based on the Kappa statistics and the correlation coefficient, indicate that SQI_W and SQI_A models provide a better evaluation of soil quality than SQI_N. The SQI_W and SQI_A models determine the quality of the soil using all the indicators and assigning higher weights to the indicators considered to be the most important in the evaluation. In contrast, the SQI_N model uses the average value of all the indicators and the lowest scoring indicator, which gives it a higher weighted value. In other words, while the SQI_W and SQI_A independently assign the score for each indicator, the SQI_N only gives preferential importance to the indicator with the lowest score.

By way of conclusion, it can be indicated that the indices that best reflect the state of soil quality are the SQI_TDS-L-W, SQI_TDS-L-A, SQI_TDS-NL-A and SQI_TDS-NL-W indices, as they are the most accurate indices and provide the most consistent results, since they consider all soil parameters. Normally, the more indicators the better the soil quality is represented, but it can happen that there is a high correlation between some of the selected properties of the soils, and this translates into a duplication of data and very laborious laboratory analysis. The evaluation method using the minimum data set avoids duplication of data and significantly reduces the labor and financial costs associated with sampling and analyzing data. As shown in the present study, the SQIMDS-L-W and SQIMDS-L-A methods provide an acceptable assessment of soil quality in large-scale studies.

Below, and as an example, the equations that define four of the quality indices that best represent the quality of soils are indicated.

As a result of applying the formula in Equation (7), the SQITDS-L-W soil quality index proposed in this work is represented in the Equation (9). The variables that make up the equation are ordered according on the weight or load value of the sixteen indicators:

SQI_TDS-L-W = (0.081 * GWC1500 kPa + 0.080 * GWC33 kPa + 0.079 * Ca + 0.076 * CEC + 0.075 * Clay + 0.073 * Mg + 0.071 * Sand + 0.068 * COLE + 0.064 * CaCO₃ + 0.062 * pH + 0.062 * EC + 0.051 * Na + 0.050 * SOC + 0.047 * Erosion factor “K” + 0.031 * BD + 0.028 * K)

(9)

If the formula in Equation (6) is applied, the SQITDS-L-A soil quality index that result can be seen in Equation (10). The variables of the equation are the sixteen scores of the indicators:

SQI_TDS-L-A = (GWC1500 kPa + GWC33 kPa + Ca + CEC + Clay + Mg + Sand + COLE + CaCO₃ + pH + EC + Na + SOC + Erosion factor “K” + BD + K)/16

(10)

The SQIMDS-L-W soil quality index is represented by Equation (11). The variables that make up the equation are ordered according to the value of the load of the six indicators:

SQI_MDS-L-W = (0.228 * GWC1500 kPa + 0.220 * Ca + 0.215 * CaCO₃ + 0.212 * Clay + 0.098 * Mg + 0.027 * SOC)

(11)

The SQIMDS-L-A soil quality index proposed in this work is given by Equation (12). The variables that make up the equation are the six indicator scores:

SQI_MDS-L-A = (GWC1500 kPa + Ca + CaCO₃ + Clay + Mg + SOC)/6

(12)

On average, the order of relative contribution of the selected indicators to the estimation of the SQI has been the following: water retention, water retention, the content of Calcium and Magnesium, the percentage of Clay, the Extractable Exchange Capacity, the content in Calcium Carbonate Equivalent, the Ph, etc. (see Equation (9)) This clearly reflected the influence of the weighting factors attributed through factor analysis.

The specific contribution of each MDS indicator to SQI shows that water retention had the highest contribution to SQI, followed by calcium content, equivalent calcium carbonate content, clay percentage, magnesium content and organic carbon content (see Equation (11)).

The validity of the method used can be confirmed with the opinion of agricultural experts in this region, who consider deep soils to be the best soils in this region (very high–high quality), which corresponds to Vertisols, Alfisols, Inceptisols and Mollisols, developed on clays, marls and deposits of glacis, the Mesozoic Covert and Sierra de Alcaraz. These are clayey soils (preferably smectite-type clays), with high organic matter content, high water retention, high CEC and COLE. The results obtained show that the consideration of the properties of the subsurface horizons of the soils and, specifically, the clay texture and water retention, are of great value for the evaluation of soil quality, considering that the climate in this region is semi-arid, and the soil moisture regime is xeric, so the soils remain in a moisture deficit for a large part of the year (between the months of May to October). The yield of rainfed crops will depend, mainly, on the moisture content stored in the subsurface horizons of the soil profile. Furthermore, the quality increases with the presence of a carbonate accumulation horizon and a pH close to neutrality, which is optimal for most crops.

According to agricultural experts, the properties that provide the highest quality to the soils of this region coincide with the indicators with the highest loads (high water retention, moderate percentage of clay, high CEC and COLE, neutral pH, moderate content of carbonates, etc., (see Equations (9) and (10)).

The soils with MODERATE quality closely correspond to the predominance of Inceptisols and Alfisols, developed on dolomites of the Sierra de Alcaraz and on deposits of the glacis that surround the Sierra del Relumbrar. Generally, with soils in which the loam–clay–sandy and clay–sandy textures predominate, in horizons A and B, respectively. In addition, they usually present clay content, water retention and a CEC with moderately high values. The pH values usually vary between 5 and 7.

The LOW quality soils correspond to thin soils (Entisols) developed on quartzites, slates and sandstones, located in the Sierra del Relumbrar, in the valleys on Triassic sandstones that surround the Sierra del Relumbrar and on dolomites in the Mesozoic Covert. of the Iberian Massif. These are soils of small thickness, formed by a single ochric epipedon, which rests on quartzites, slates, dolomites or sandstones. They are noncalcareous soils; the texture is usually loam–sandy or sandy–loam. The pH values vary between 4 and 6. They are very poor soils in nutritive elements for plants, with very low water retention; the degree of saturation is very close to 50% and the COLE has very low or insignificant values. These soils, traditionally, have been lent for land use dedicated to extensive livestock (grasslands).

3.5. Soil Quality Assessment Based on Land Use and the Different Parent Materials

The resulting values corresponding to the soil quality indices TDS-L-W and TDS-L-A of the three different types of land use (woodland, grassland and cropland) were 0.45, 0.45 and 0.52, respectively. The values of the TDS-NL-W and TDS-NL-A indices were 0.41, 0.44 and 0.50, respectively. The values of the indices calculated by the Nemoro method (TDS-L-N and TDS-NL-N) were 0.29, 0.30 and 0.35, respectively.

In summary, the TDS indices calculated using a linear and nonlinear equation and with the additive and weighted additive methods give higher values than the same indices, but calculated by the Nemoro method.

Similar results were obtained with the MDS indices calculated by means of a linear and nonlinear equation; with the additive and weighted additive methods they gave higher values (between 0.35 and 0.54) than the same indices, but calculated by the Nemoro method (between 0.24 and 0.34).

The mean soil quality values under crops were all significantly higher than those of the other land use. Due to the higher values of F and p (Figure 8), the values of the indices calculated with a linear equation were more sensitive than those calculated with a nonlinear equation.

The mean values for the twelve quality indices used were significantly higher for farmland than for grassland and forest soils. This fact suggests that there is a correlation (highly significant (p < 0.01 **) or probably significant (p < 0.05 *)) between the quality index and the land uses. The mean values of soil quality under forest use have been found to be equal to or slightly lower than that of grassland use.

Traditionally, the farmers of this region in this region have selected soils with the best edaphic properties for their agricultural use, and although, over time, conventional agriculture causes the degradation of some soil properties, agricultural soils and anthropogenic soils have maintained higher-quality levels compared to natural soils. In contrast, soils located in areas with a steep slope or shallow soils, with abundant rocky or stony, sandy texture, with low CEC, etc., are those that have been used solely and exclusively for grassland and/or forest use (holm oak cleared up).

In addition, in Figure 8 it can be seen that the lowest average values of the quality indices are those obtained by the Nemoro method, both using linear and nonlinear scoring.

The resulting values corresponding to the quality indices TDS-LW and TDS-LA, of the four different types of source rocks of the soils (quartzite, sandstone, limestone and clay–loam) were 0.41, 0.45, 0.49 and 0.55, respectively. The values of the TDS-NL-W and TDS-NL-A indices were 0.36, 0.41, 0.48 and 0.52, respectively. The values of the indices calculated by the Nemoro method (TDS-L-N and TDS-NL-N) were 0.25, 0.29, 0.33 and 0.36, respectively.

Similar results were obtained with the MDS indices calculated by means of a linear and nonlinear equation; with the additive and weighted additive methods they gave higher values (between 0.23 and 0.58) than the same indices, but calculated by the Nemoro method (between 0.17 and 0.35).

The behavior of the soil quality indices in relation to the bedrock has been studied (Figure 9); the mean values for the twelve quality indices used were significantly higher for soils formed on clays and/or loams, while soils formed on quartzites and/or blackboards presented the lowest SQI. The soils developed on limestones and/or dolomites presented a quality index that we could consider as moderate–high and moderate–low quality values were obtained on soils formed on sandstones. There is a highly significant correlation (p < 0.001 ***) between all the quality indices and the parent rocks of the soils.

Again, it can be seen in Figure 9 that the lowest average values of the quality indices are those obtained by the Nemoro method. In this sense, it can be observed in Figure 8 and Figure 9 that the highest values of the quality index are those obtained through linear and nonlinear equations.

The quality indices obtained by the MDS method evaluate the soils of the studied region through a selection of only six indicators, which represent the highest loads obtained through factor analysis (FA). Well, the order of relative contribution of these indicators to the estimation of the soil quality index is as follows: 22.8% for water retention, 22% for extractable Ca, 21.5% for CaCO₃ equivalent, 21.2% for the clay content, 9.8% for the extractable Mg and 2.7% for the OC stock.

The score values of the water retention indicator at field capacity (GWC_1500kPa) varied from 0.41, 0.50 and 0.58 (linear) and 0.36, 0.45 and 0.52 (nonlinear) for grasslands, forest soils and farmland, respectively (Figure 10A,B). These indicators score values (GWC_1500kPa), both linear and nonlinear, were highly significant (p < 0.01).

A similar pattern was observed with the indicator score values (Clay), since the highest values are those obtained using a linear equation and correspond to croplands. On the contrary, the lowest correspond to grasslands

The score values of the indicators (CaCO₃) and (Ca) have a totally different pattern from the previously studied indicators, since the highest values are those obtained using a nonlinear equation corresponding to croplands and the lowest values to forest soils. The score values of the indicators (OC Stock) and (Mg) also present as higher values those obtained using a nonlinear equation, but the highest values usually correspond to forest soils and values of pastures and farmland. They are very similar and even the same. These score values of the indicators (OC Stock) and (Mg), both linear and nonlinear, were significantly lower (p < 0.05).

The score values of the water retention indicator at field capacity (GWC_1500kPa) varied from 0.39, 0.42, 0.50 and 0.69 (linear) and 0.33, 0.36, 0.45 and 0.63 (nonlinear) for soils developed on sandstones, quartzites, limestone–dolomites and clays–marls, respectively (Figure 10C,D). These indicators score values (GWC_1500kPa), both linear and nonlinear, were highly significant (p < 0.001).

A similar pattern has been observed with the indicator score values (Clay), since the highest values (0.61, 0.66, 0.71 and 0.82) are obtained using a linear equation. Additionally, these higher values correspond to the soils developed on clay–loams and the lowest to those formed on sandstones.

The linear score values of the indicator (OC Stock) indicate that the soils with the highest content of organic matter are the soils developed on clay–loams and those with the least content are those formed from sandstones. On the contrary, the nonlinear score values offer us a different pattern, since the soils with the highest content of organic matter are the soils developed on sandstones and those with the least content are those formed from clay–marl and limestone–dolomites.

The score values of the indicators (CaCO₃) varied from 0.00, 0.17, 0.32 and 0.39 (linear) and 0.00, 0.30, 0.55 and 0.58 (nonlinear), for soils developed on quartzites, sandstones, clays–marl and limestone–dolomites, respectively.

A pattern similar to the previous one has been observed with the indicator score values (Ca).

The highest indicator scores values (Mg) have corresponded to soils developed on clay–loams and the lowest values to those formed from sandstones.

Figure 11A shows the order of specific contribution of the six selected indicators to the estimation of the soil quality index (expressed as averages), according to the MDS method, for different land uses, using linear scoring functions. Between 60 and 70, 60% of the contribution of the indicators to the quality of forest soils and croplands is constituted, in decreasing order by (Clay), (GWC_1500kPa) and (Ca); in soils under grassland, the highest contribution (80%) is constituted by (Clay), (OC Stock), (GWC_1500kPa) and (Ca). The indicators (Clay), (GWC_1500kPa) and (Ca) can be considered as having the greatest influence on soil quality. In contrast, (Mg) and (CaCO₃) are the least influential.

However, using non-linear scoring functions (Figure 11B), 80% of the contribution of the indicators to the quality of forest soils is made up, in decreasing order, by (OC Stock) (Mg), (GWC_1500kPa) and (Clay); in grassland soils the largest contribution (70%) is made by (OC Stock), (Clay), (Ca) and (GWC_1500kPa) and in cropland by (Ca), (CaCO₃), (OC Stock) and (GWC_1500kPa). In this case, in the non-linear scoring functions, the indicators with the highest influence on soil quality are (OC Stock), (Ca), (GWC_1500kPa) and (Clay) and the indicator with the lowest contribution to soil quality is (Mg).

Figure 11C shows the order of specific contribution of the selected indicators (expressed in averages) to the estimation of the soil quality index, for soils developed on different mother rocks, using linear scoring functions. Eighty percent of the contribution of the indicators to the quality of the soils developed on quartzites is constituted by (Clay), (GWC_1500kPa) and (OC Stock), in soils on sandstones and limestone–dolomites and clays-marls by (Clay), (Ca), (GWC_1500kPa) and (OC Stock). (Clay), (GWC1500 kPa), (Ca) and (OC Stock) can be considered as the indicators with the greatest influence on soil quality, and (Mg) and (CaCO₃) as those with the least contribution.

However, by using nonlinear scoring functions (Figure 11D), 90% of the contribution of the indicators to the quality of the soils developed on quartzites is constituted by (OC Stock), (Clay), (GWC_1500kPa) and (Mg); in soils on sandstones the greatest contribution is formed by (OC Stock), (Ca), (Clay) and (GWC_1500kPa); in soils formed from limestone–dolomites by (Ca), (CaCO₃) and (OC Stock) and in soils on clay–loams by (Ca), (GWC_1500kPa), (CaCO₃) and (OC Stock). In this case, the indicators with the greatest influence on soil quality are (OC Stock), (Clay), (Ca), (GWC_1500kPa), and the indicator with the least contribution to soil quality is (Mg).

In summary, the indicators with the greatest influence on soil quality, obtained by the MDS method, for different land uses and those developed on different rocks, using linear scoring functions, are the following: (Clay), (GWC_1500kPa) and (Ca). These data, in turn, coincide with the opinion of the agricultural experts of this region (who affirm that the best soils in this region are deep soils, with a clayey texture, with high water retention and a pH close to neutrality). However, the indicators with the greatest influence on soil quality, using nonlinear scoring functions, are: (OC Stock), (Clay), (Ca) and (GWC_1500kPa). In other words, the most important indicator is the OC stock, which is not logical in the case of a region in which the soils have a very low average soil OC content (0.86%).

4. Conclusions

With the research carried out in this study, it has been possible to identify the indicators that have a greater implication or load in the quality index of the soils of Campo de Montiel-Sierra de Alcaraz, which are, in descending order, the following: water retention, extractable calcium, extractable exchange capacity, clay content, extractable magnesium, sand content, COLE, calcium carbonate equivalent content, pH and electrical conductivity. On the contrary, the indicators that have a lower implication or load in the quality index are: extractable potassium, apparent density, K factor (erosion), organic carbon stock and extractable sodium.

The first objective of this work consisted of evaluating the quality of the soil, focusing on the support function of natural ecosystems and the production of food and fibers in the studied region, and as a conclusion it has been obtained that the indices that best reflect the state of the quality of the soil, in studies at a detailed or semi-detailed scale, are the SQI_TDS-L-W, SQI_TDS-L-A, SQI_TDS-NL-W and SQI_TDS-NL-A indices (indices made with all the indicators). In addition, the SQI_MDS-L-W, SQI_MDS-L-A, SQI_MDS-NL-W and SQI_MDS-NL-A indices provide an acceptable assessment of soil quality in large-scale studies (indices performed with the minimum data set). The results obtained in this study indicate that the linear and nonlinear scoring equations (L and NL) used to estimate soil quality seem to work equally well, both with the Weighted Additive (W) and Additive (A) methods. However, the same is not the case with the Nemoro (N) model. That is, the results obtained indicate that the SQI_A and SQI_w indices could provide a better estimate compared to SQI_N.

Another of the objectives of this work was to make maps with the quality indices obtained, in order to identify the worst and best-quality areas, and as a conclusion it wasfound that the maps showing the spatial distribution of the different grades of soil quality show that the moderate and high-quality soil areas were dominant and represented around 70% of the area in the studied region. Very high and very low-quality areas are very limited. The surface percentages of the different soil quality grades obtained, in the TDS and MDS models, are similar with the use of the Weighted Additive (W) and Additive (A) methods; however, with the Nemoro (N) method the results are inferior (it can be said that they are underestimated).

In order for the evaluation of soil quality to provide us with the most complete information possible, this evaluation of the quality indices must be carried out considering the properties of all the horizons of the soil control section (between 0 and 100 cm deep).

Another objective of this work was to establish a relationship between the quality indices obtained with the uses of the land and with the mother rocks of the soils, and as a conclusion it was found that the quality indices were higher on croplands than on grasslands and forest soils. In turn, the highest quality indices were obtained in soils developed on marls and clays, while soils formed on quartzites and slate present the lowest indices.

The A horizons of forest and grassland soils have an average content of 3.5% of organic carbon, while the surface horizons of cultivated soils have an average content of 0.8%. This difference shows the rapid degradation of the organic matter in the surface horizons due to the cultivation effect. Under climate change conditions, Spanish agricultural soils could act as potential carbon sinks considering that the potential capacity for additional OC sequestration is higher in cropland soils than in grassland soils due to their deeper soil profiles and lower OC saturation. For this, it is necessary to promote agricultural techniques that favor the conservation and increase of carbon, such as conservation tillage, addition of exogenous organic matter and cover crops and fallows with vegetation.

The intensive and continuous cultivation, the elimination of residues, the low application of organic fertilizers—which leads to the decrease in organic matter—and the degradation (erosion) are the main causes of the reduction of the productivity of the soil in the systems of rainfed cultivation in this semi-arid region of the Iberian Peninsula.

The degradation of soil quality with excessive or conventional tillage is particularly important in agroecosystems located in arid and semi-arid environments. Agricultural management practices using reduced tillage or conservation systems should be proposed as a means to maintain or even improve organic carbon storage and soil quality in arid and semi-arid regions [55,63]. This reduced tillage causes less physical disturbance in the soil and improves the quality of the soil by improving the structure of the soil, by increasing the content of organic carbon and the microbial community of the soil, which is important for the transformation and mineralization of organic compounds and nutrients in the soil ecosystem [64,65].

The yield of rainfed cereals (wheat and barley) in this Mediterranean region depends mainly on the available moisture stored in the soil profile during the rainy season. In this semi-arid ecosystem, annual evapotranspiration is much higher than precipitation. Therefore, soils with higher clay content, in turn, have higher water retention, higher CEC, etc., are considered the highest quality in the region and are of great importance for maintaining agricultural productivity.

It has been widely recognized that soil carbon sequestration can be of great importance as a mitigation and adaptation measure to climate change. However, it is often forgotten that organic carbon plays an equally important role in ensuring food security. Climate change is likely to have a strong impact on agriculture, posing a major threat to food security.

The methodology presented in this study, carried out in a semi-arid ecosystem, could be applied in other regions with different ecosystems, although surely the indicators with a greater implication or load in the soil quality indices are different from those obtained here.