Spatial Variations of Vegetation Index from Remote Sensing Linked to Soil Colloidal Status

Abstract

Recent decades have seen a progressive degradation of soils owing to an intensification of farming practices (weeding and high trafficking), increasing use of pesticides and fertilizers, mainly nitrogen, resulting in a steady decline in soil organic matter, a key component to maintain soil fertility. The work has coupled the normalized difference vegetation index (NDVI) of wheat cultivation in Central Italy to soil properties where the wheat was grown to identify the properties linked to within-field variability in productivity. NDVI was assessed through Copernicus Sentinel-2 (S-2) data during the wheat anthesis phase. The main outcome showed a significant correlation of NDVI variability to soil colloidal status and to the relative quantity in the exchange complex of the Ca${}^{2+}$ ions. No relationship emerged between NDVI and soil macronutrients (nitrogen, phosphorus, and potassium) concentration. The work suggested that such elements (nitrogen, especially) should not be provided solely considering the vegetation index spatial variations. Rational and sustainable management of soil fertility requires the integration of the NDVI data with the whole complex of soil physical/chemical status. In this way, the identification of the real key factors of fertility will avoid the negative impact of overfertilization. As an example, a fertilization plan was simulated for the sunflower–wheat sequence. The results showed that in the study area additional supplies of N and K would be unnecessary.

1. Introduction

Approximately 81% of the organic carbon resources that are actively involved in the global carbon cycle are stored in soils [1]. Soil organic matter (SOM) represents one of the largest reservoirs of carbon on the global scale; its quantity and quality are important in the management of soil fertility, nutrient supply, and carbon dynamics [2]. Preserving and/or increasing the SOM pool ensures favorable nutrient conditions for field crops and, in turn, contributes to securing food security [3,4]. It has been remarked that soil degradation should be recognized, alongside climate change, as one of the most pressing problems facing humanity particularly in arid and semi-arid regions where salt-induced soil degradation coupled to intensive farming is a major cause of soil organic carbon (SOC) loss [5], the main component of SOM.

Although poor irrigation tilling practices coupled to the extensive use of pesticides and fertilizers represent a major cause of soil organic matter loss, degradation of the fertile layer of the soil, and its desertification [6,7], a few farming practices can be employed to improve soil functional qualities and increase SOC, including optimal fertilization, crop–grassland rotation, and amendments [5].

In the Mediterranean area, a typical agriculture soil is characterized by a high content of clay and active limestone; declining SOM; and by a high level of degradation due to intensive use, frequent mineral fertilization, weeding, no addition of organic matter, and intense traffic of heavy-weight machinery [6].

Available multi-layered studies associate vegetation indexes from remote sensing to nitrogen fertilization [8], to the content [9] or to the supply of N, P, and K [10], to soil texture [11], or soil type [12]. Less frequent (or absent) is the association between remote data and specific soil chemical and physical features.

Remote sensing technologies have proven to be viable for assessing and mapping the availability of agro-forestry resources such as productivity, residual biomass, and crop yields [13] as well as to optimize the interaction crop–soil environment with the final objective of increasing the sustainability of yield by means of more judicious input management [14,15]. Focused mainly on site-specific fertilization [16,17], remote sensing application has progressively been extended to the analysis of all the sources of intra-field variation [18,19] such as yield, soil, crop, anomalous factor (i.e., weeds or pathogens), and management variability (from tillage to the application of fertilizers or pesticides).

Depending on the properties of the instrument and platform, remote sensing data are available from coarse (more than one kilometer) to fine (sub-meter) spatial resolution, and at variable temporal resolution, daily to monthly [20]. The main remote and proximal sensing technologies employed in precision agriculture include satellite platforms, drones, and sensors installed on tractors [15]. Earth observation satellites record and collect spatial information regularly, with wide coverage and low cost, and therefore represent an advantageous tool for the detection of natural and agricultural resources over the last decades [21]. Sentinels 1 and 2 or lower resolution satellite missions have been used to create dynamic cropland masks [22]; perform crop type mapping [23]; estimate soil moisture [24,25]; monitor rice production [26]; estimate plant parameters such as leaf area index [27,28]; or aid in within-field decisions in precision agriculture settings for crop yield mapping, fertilizer use, and minimizing nitrogen loss to water [8,29,30].

We hypothesized that the mineral status of the soil layer where the great part of wheat roots are distributed, due to its influence on plant nutritional status [31,32], should be picked up by remotely sensed normalized difference vegetation index (NDVI) layer.

An exploratory investigation was applied to infer a fundamentally general relationship among exploratory (soil mineral status) and response variables (NDVI) from remote sensing in a typical Mediterranean field. A wheat cultivation was set up in late 2018 to evaluate the contribution of NDVI, calibrated by soil chemical properties monitoring (Figure 1). The study attempted to answer two key questions: (1) which soil nutritional elements and physical features account for the within-field spatial variation of crop productivity, and (2) whether and how the NDVI can be used as a proxy of soil mineral status.

2. Methods

2.1. Study Site

A 5.0 ha experimental field crop, spatially located in the northeastern outskirts of Rome, central Italy (Figure 2, Latitude 42.103° N, 12.628° E), was set up in 2018. The field is property of the Research Center for Engineering and Agro-Food Processing (CREA) to pursue academic investigations.

Geologic formation outcrops throughout the area consist of volcanic tuff effusive types of lower and middle Pleistocene. These formations are connected to the intense volcanic activity of the northern Lazio region at the time, particularly associated with the calderas activity of the area. The predominant rock types are pyroclastic launch products, mainly composed of loose sand-lapilli levels and sometimes with the presence of cineritic levels more or less cemented. The pyroclastics are leucititic type, which is of significant importance in the process of alteration, show an intense activity of quaternary hydrothermal type, which led to the formation of Analcime (natural zeolite). This confers special properties to the soil in terms of water retention and nutrients release (less water and nutrients available in the real conditions with respect to the analytical data, as captured by the zeolite lattice structure). The southern part of the experimental farm of CREA also includes a portion of an alluvial valley of a secondary stream of the Tevere River. The alluvio-colluvial deposits here consist of fine sandy loam and fine sediment resulting from erosion and reworking of the deposits and soils of the slopes [33]. Soils are of volcanic origin and are classified as Typic Argixeroll [34], soil profile was described in [35].

From the climatic and pedoclimatic point of view, according to the long-term data (30 years) of the ISIS 1.0 database, the average annual air temperature is 13.7 °C, the average annual rainfall is 890 mm, equivalent to an aridity index calculated with the De Martonne equation of 37.5 class (35–40) moderately humid. The thermal regime of the soil is thermal (15–22) with an average soil temperature of 16.3 °C at 0.50 m depth. The water regime is xeric (80–115 days), with 88 cumulated dry days per year. Climatic data were collected by the monitoring station adjacent to the study area [33].

The soil being tested was characterized through a complete physical-chemical analysis, in order to evaluate the characteristics of the mineralogical components and the relationships between them; the reactions and the electrical conductivity that influence the bioavailability of many nutrients; the level of chemical fertility, through the determination of the nutrient content and the colloidal capacity of the soil, which indicates the state of aggregation of the particles; the drainage ratio; and the buffering power of the soil. The analysis involved the collection within the area of twenty random sampling points considered as replicates. The sampling points were randomly drawn in a downsized study area to account for a 10 m border effect (Figure 2).

A 2 kg soil sample from the 0–20 cm layer in each sampling point was collected. At wheat heading 50% of roots are localized in the 0–20 cm soil layer while a further 10% being found in the 20–40 cm layer [36]. The figures were confirmed by the authors of [37] by modeling root distributions of eleven temperate crops: at least half of the root biomass could be found in the upper 20 cm of soil, 61–68% of wheat roots are found in the 0–30 cm soil layer.

Soil characterization was carried out according to the Italian official method of analysis [38] by a UNI CEI EN ISO/IEC 17025:2005-certified laboratory. Particle size distribution was determined by gravimetric method; pH in H${}_2$O with a potentiometer; total organic carbon through Walkley and Black’s method; total nitrogen with Kjeldahl’s method; available phosphorus by means of Olsen’s procedure; cation exchange capacity and exchangeable cations measured in the extracted soil solution (ammonium acetate) by using the atomic absorption spectroscopy (AAS, model AA240FS, Varian, Crawley—UK); assimilable metals by extraction with diethylenetriaminepentaacetic acid (DTPA); and spectrometric analysis with a inductive coupled plasma spectrophotometer model iCAP Pro (ICP-OES, Thermo Fisher Scientific, Waltham—USA). Soil analyses included macro- and micronutrients as well as soil particle size distribution (

Table 1. Soil analyses and associated units of measure carried out on the 20 sampling set points. Percentages are mass over mass (m/m). Abbreviations: Tot., Total; Act., Active; Av., Available; Ass., Assimilable; Sol., soluble; Exc., exchangeable; CEC, Cation Exchange Capacity; EC, Electrical Conductivity; CI, Colloids Index; SOM, Soil Organic Matter.

Physical Properties	Chemical Properties	Derived Chemical Properties
Sand [%]	Tot. and Act. limestone [%]	Magnesium/potassium ratio
Silt [%]	Tot. organic carbon [%]	Carbon/nitrogen ratio
Clay [%]	Tot. nitrogen [%]	SOM [%]
Particles size distribution [%]	Av. phosphorus [mg kg${}^{-1}$]	CI [%]
pH	Ass. iron [mg kg${}^{-1}$]	Exc. potassium [% of CEC]
EC [mS cm${}^{-1}$]	Ass. manganese [mg kg${}^{-1}$]	Exc. magnesium [% of CEC]
	Ass. copper [mg kg${}^{-1}$]	Exc. sodium [% of CEC]
	Ass. zinc [mg kg${}^{-1}$]	Exc. calcium [% of CEC]
	Sol. boron [mg kg${}^{-1}$]
	Sol. cobalt [mg kg${}^{-1}$]
	Exc. calcium [mg kg${}^{-1}$]
	Exc. magnesium [mg kg${}^{-1}$]
	Exc. potassium [mg kg${}^{-1}$]
	Exc. sodium [mg kg${}^{-1}$]
	CEC [meq 100 g${}^{-1}$]

).

According to the US classification standards, soil particles are divided into three grades: clay particles <2 $\mathsf{\mu }$m, silt particles ≥2 $\mathsf{\mu }$m <50 $\mathsf{\mu }$m, and sand ≥50 $\mathsf{\mu }$m <2000 $\mathsf{\mu }$m.

Among the derived chemical properties, exchangeable potassium, magnesium, sodium, and calcium are defined in % of the CEC; colloids index is defined according to [39]

$$CI=10·SO{M}_{\%}+Cla{y}_{\%}$$

(1)

2.2. Vegetation Model

At the end of 2018, the study area was sown with durum wheat. The cultivation was conducted following the common farming practices of the area. A linear model was set up to investigate the effect of chemical and physical properties of the soil and soil/vegetation moisture on the Normalized Difference Vegetation Index (NDVI) of the wheat crop growing during winter and spring of 2019, on the same 20 sampling points identified for soil analysis. NDVI is directly related to the photosynthetic capacity and therefore to the energy absorption of plant canopies [40,41], thus proving to be an excellent predictor of productivity and yield [30].

NDVI was estimated for the study area on all available passes of the satellites of the Copernicus Sentinel-2 (S-2) mission in 2019 as

$${NDVI}_{t,x,y}=\frac{{\rho }_{842,t,x,y}-{\rho }_{490,t,x,y}}{{\rho }_{842,t,x,y}+{\rho }_{490,t,x,y}}$$

(2)

where $NDV{I}_{t,x,v}$ is NDVI at time t and spatial coordinates x, y; ${\rho }_{490,t,x,y}$ and ${\rho }_{842,t,x,y}$ are the spectral reflectances of the central wavelengths of the near-infrared and red bands of S-2 recorded at time t and at x, y coordinates. These spectral reflectances are themselves ratios of the reflected over the incoming radiation in each spectral band.

The S-2 satellites aim at providing multispectral data with a 5-day revisit frequency and 10 meters spatial resolution [42]. The medium-to-high spatial resolution granted the independence assumption of the sampling points locations (i.e., a one-to-one correspondence links the soil sampling point set and the S-2 pixel set). Cloud and cirrus formations were detected and removed through the quality assurance metadata provided and the resulting pixels masked from NDVI calculation.

The NDVI profile of the study area (Figure 3) helped in tracing the timing of phenology of the crop. At the latitudes of the study, wheat sowing takes place between the end of October and the beginning of November. Field observation [43], phenological model [44,45], or analysis of vegetation index [46] confirmed that in Mediterranean environments, the anthesis occurs between the end of spring frosts and the beginning of the summer drought, corresponding to the end of April–first half of May. The passage from anthesis to maturity is a crucial phase shift because the photosynthates accumulated in the photosynthetic organs (source) relocate towards the ear (sink) for grain filling [47,48]. NDVI was further calculated for the 20 sampling points set, on the peak season day image, this $NDV{I}_{x,y}$ variable was used as a model predictor. The estimation of the NDVI profiles was performed in Google Earth Engine [49].

Soil chemical and physical features on the 20 sampling points set along with vegetation and soil moisture were included in the model as potential explanatory variables. Moisture was proxied by extracting the C-band Synthetic Aperture Radar (SAR) Ground Range Detected (GRD) single bands on the sampling points for the Sentinel-1 image available on 15 May 2019. Vertical–Vertical (VV) and Vertical–Horizontal (VH) bands report the portion of the outgoing radar signal that the target redirects directly back towards the radar antenna; in VV mode, the microwaves of the electric field are oriented in the vertical plane for both signal transmission and reception whereas in VH mode the backscatter signal is received in the horizontal plane.

The image was preprocessed on Google Earth Engine (apply orbit file, GRD border noise removal, thermal noise removal, radiometric calibration, and terrain correction). The process to define the optimal minimal model of NDVI prediction from the starting comprehensive model formed by all explanatory variables included four steps. First, all variables were standardized (i.e., centered on their mean and scaled by their standard deviation) to ease interpretation of model coefficients and avoid issues due to multicollinearity of explanatory variables. In these standardized models, a unit increase in an explanatory variable is equal to its standard deviation, and it affects the predictive variable by a unit of its standard deviation. Intercept term was dropped from the standardized models of the successive steps.

Second step concerned feature selection of explanatory variables. Firstly highly correlated variables were removed. If two variables had a high correlation (pairwise absolute correlation cutoff: 0.95), the variable with the largest mean absolute correlation was removed. Highly correlated removed variables included carbon/nitrogen ratio, organic carbon, SOM, and exchangeable calcium (mg kg${}^{-1}$) (Figure 4). Second, important variables were selected by fitting an unsupervised random forest classification model over different tuning parameters and filtering out the least significant variables based on the importance measure, on a percentage scale (cutoff: 40%). Among the remaining explanatory variables, Colloids Index (CI), CEC, exchangeable Calcium (% of CEC), and total nitrogen were selected in order of decreasing importance. This intermediate linear model explained 67% of variance in NDVI (Adjusted R-squared: 0.6782); residual standard error was 0.55 on 16 degrees of freedom.

The third step involved stepwise model selection, based on Akaike’s An Information Criterion (AIC) value, of a set of linear models fitted using generalized least squares. Each variable was considered for subtraction from the set of explanatory variables based on AIC value.

Heteroskedasticity and spatial correlation were accounted for in the fourth step. A slightly increasing linear relationship in residuals vs. NDVI values for the explanatory variables was accounted for by weighting observations by selecting variance functions that minimized AIC while being not significantly different from the optimal model. Variance functions chosen were

• an exponential function for CI (where denotes the variance function evaluated at CI and t is the variance function coefficient, t = −0.39);
• a power function for exchangeable Ca (t = −0.32).

Similarly, spatial autocorrelation was accounted for by evaluating the better performing correlation structure (longitude + latitude) in terms of AIC that resulted the spherical spatial correlation (where d is the range and n is the nudge, d = 91.2, n = 0.004). Model estimation was performed in R 3.6.3 [50], stepwise procedure by package MASS 7.3–51.5 [51], GLS modeling by package nlme 3.1–144 [52], variable importance by packages caret 6.0–85 [53], and randomForest 4.6–14 [54]; general data table management was performed by package data.table 1.12.8 [55].

2.3. Fertilization Plans

For applying the procedure to a case study, nitrogen, phosphorus, and potassium fertilization plans were computed from soil chemical and physical analyses of the 20 sampling points set, following the indications formulated in the regulation drawn up by the Lazio Region [56]. Soil nutrient balances taking into account inputs and losses of N, P, and K were computed to satisfy the nutrient demands of sunflower and wheat. As an example, nitrogen balance (a dynamic element considered fundamental for productivity) included seven components: N crop demand for sunflower and for wheat, availability of N for the crop, N leakage caused by rainfall, N leakage due to immobilization processes, residual N supply from previous crop, residual N supply from previous organic fertilizations, and supply of N from natural and anthropic sources. Expected yields were estimated from 2019 average yields in the province of Rome (http://dati.istat.it/): 1330 kg/ha for sunflower and 3000 kg/ha for wheat. The fertilization plan was estimated by fertplan 0.1 [57], an R package specifically developed; spatialization of the fertilization plan from the sampling set was performed by ordinary kriging in R 3.6.3 by package gstat 2.0–4 [58].

3. Results

All the analyzed soil samples fall within the clayey loam or silty clayey loam USDA classification (

Table 2. Descriptive statistics for remote sensing and soil physical properties on sampling point set. Abbreviations: AVG: average, MIN: minimum value, MAX: maximum value, STD: standard deviation, CV: coefficient of variation (i.e., $STD/AVG·100$); NDVI: Normalized Difference Vegetation Index. See Table 1 for soil physical abbreviations and units of measures and Figure 4 for remote sensing abbreviations.

	Remote Sensing Properties			Soil Physical Properties
	NDVI	VH	VV	Moisture	Sand	Silt	Clay	pH	EC
	$\times {\mathbf{10}}^{\mathbf{4}}$	Backscatter		%	%	%	%		mS cm${}^{-\mathbf{1}}$
AVG	9288	−21.4	−15.08	20.5	25	38	37	7.9	0.47
MIN	9124	−22.9	−16.72	18.2	16	21	34	7.4	0.32
MAX	9404	−19.6	−13.16	23.8	41	48	40	8.0	0.65
STD	76.27	0.98	0.89	1.29	7.0	7.6	2.2	0.10	0.06
CV	0.82%	−4.50%	−5.90%	6.28%	28%	20%	5.9%	1.7%	13%

; Figure 5) and have a sub-alkaline reaction. Workability is difficult with a tendency to retain too much water, often resulting in stagnant water after heavy rainfall. Water stagnation can be a serious problem for most crops. It often causes stunted growth and rot or other diseases. Therefore, this type of soil requires drainage of surface waters.

Total limestone is high, with a high content of active limestone (

Table 3. Descriptive statistics for soil chemical properties on sampling point set. See Table 1 and Table 2 for abbreviations and units of measures.

		AVG	MIN	MAX	STD	CV
Tot. lime	[%]	9.69	0.6	17.4	4.01	41.40%
Act. lime	[%]	4.30	0.0	6.5	1.55	36.00%
SOM	[%]	2.63	1.97	3.56	0.352	13.40%
Org. C	[%]	1.52	1.14	2.06	0.204	13.40%
Tot. N	[%]	0.157	0.122	0.205	0.019	11.80%
Ass. P	[mg kg${}^{-1}$]	14.1	11.0	21.0	2.57	18.20%
Ass. Fe	[mg kg${}^{-1}$]	23.4	16.4	32.0	3.99	17.10%
Ass. Mn.	[mg kg${}^{-1}$]	36.3	22.8	70.4	12	33.00%
Ass. Cu	[mg kg${}^{-1}$]	4.91	4.0	6.40	0.676	13.80%
Ass. Zn	[mg kg${}^{-1}$]	1.31	0.60	3.40	0.626	48.00%
Exc. Ca	[mg kg${}^{-1}$]	5852	4940	6640	390	6.67%
Exc. Mg	[mg kg${}^{-1}$]	149.5	124	218	20.2	13.50%
Exc. K	[mg kg${}^{-1}$]	365.3	285	492	58.7	16.10%
Exc. Na	[mg kg${}^{-1}$]	87.95	66	124	14.2	16.20%
Exc. Ca	[% of CEC]	91.71	89.17	92.93	0.909	0.99%
Exc. Mg	[% of CEC]	3.875	3.254	6.141	0.588	15.20%
Exc. K	[% of CEC]	2.953	2.174	4.211	0.549	18.60%
Exc. Na	[% of CEC]	1.201	0.894	1.637	0.173	14.40%
Sol. B	[mg kg${}^{-1}$]	0.872	0.5	1.36	0.21	24.00%
Sol. Co	[mg kg${}^{-1}$]	0.014	0.01	0.02	0.005	36.20%
CEC	[meq 100 g${}^{-1}$]	31.8	27.3	35.8	1.95	6.11%
Mg/K		1.37	0.9	2	0.27	19.70%
C/N		9.71	9.34	10	0.153	1.58%
CI	[%]	63.2	57	71.6	3.7	5.85%

). The cation exchange capacity is high (>30 meq 100 g${}^{-1}$), as the basic cations saturation rate and colloids index. Calcium is the most present exchangeable cation, and the activity of limestone causes this metal to completely saturate the exchange complex. The content of organic matter and total nitrogen can be classified as medium; given the sub-alkaline nature of the soil, potassium is also well represented. The available phosphorus is low, highlighting a high degree of immobilization of the element due to the excess of Ca and the presence of active limestone. Metallic microelements are represented in good quantities, in particular as regards Fe and Mn. Among the macronutrients, N and K concentrations are widely above the levels of sufficiency for an average demanding crop, so their supply is not required. Phosphorus, on the other hand, is affected by the high degree of immobilization of the soil, so it must necessarily be added with specific fertilizations. In addition to water, soil also retains nutrients, greatly increasing its chemical fertility. The presence of calcium carbonates derived from the degradation of the original or secondary minerals, associated with mineral or organic colloids, contribute to the formation of a stable structure. On the other side, the release of sodium from sodium salts, represents a highly destructuring factor.

These soils, very common in the Mediterranean area, are generally used for the cultivation of cereals, oil, or industrial crops, resulting in a massive use of fertilizers and pesticides to prevent a decrease of their fertility. Their tendency to lose their structure can be contrasted with the addition of organic matter to counteract the excessive presence of clay and silt.

The NDVI from Sentinel-2 optical bands is very close to the theoretical maximum ($1\times {10}^4$) and with fairly low variability among the sampling set (Coefficient of Variation 1%, Table 2), whereas microwave bands from remote sensing exhibit higher spatial diversity.

Colloids Index (CI) and exchangeable Calcium (% of CEC) were the explanatory variables selected in the optimal vegetation model after pruning of all the other explanatory variables:

$$NDVI=f({\beta }_1·CI+{\beta }_2·Ca)$$

(3)

Standardized beta coefficients along with standard deviation of the explanatory variables are given in

Table 4. Vegetation model results. $\beta$ coefficients and their standards errors in standardized form together with their significance values and cross-correlation. $\beta$ coefficients are converted to unstandardized form by multiplying their variable standard deviations. Colloid Index is defined in Equation (1). Exc. Ca. is defined as % of Cation Exchange Capacity (Table 1). Abbreviations: Expl., Exploratory; Var., Variable; std. err., standard error; std. dev., standard deviation.

Expl. Var.	Standardized $\mathit{\beta } \pm$ Std. Err.	t-Value	p-Value	Correlation	Var. Std. Dev.	$\mathit{\beta } \pm$ Std. Err.
CI	$0.66\pm 0.04$	15.91	0		3.7	$2.5\pm 0.15$
Exc. Ca	$0.29\pm 0.05$	5.4	0	−0.08	0.91	$0.27\pm 0.01$

. Residual standard error of the optimal model decreased to 0.31 on 18 degrees of freedom while ANOVA confirmed it to be not significantly different from the intermediate model, built after variable selection. Neither of the three soil macronutrients (N, P, and K) nor carbon and any of the micronutrients were able to explain variation of NDVI in the field, during the key phenological phase of anthesis-start of grain filling of the wheat grown in 2019.

Nevertheless, among the sampling point set, although the NDVI range is very narrow (Figure 6), its variability is largely related to soil colloidal status (CI) and, to a lesser extent, to relative quantity in the exchange complex of the Ca${}^{2+}$ ions. Colloids index, in turn, heavily depends on SOM variations (×10) coupled to clay soil quantity [39] so that a limited increase in SOM can greatly improve soil colloid status leading to an increase in NDVI.

The spatialized fertilization plans and the relative doses, calculated in accordance with the regional guidelines, showed a different pattern of within-field variability for N, P, and K. The maps confirmed as N and K concentrations were above the demand for both sunflower and wheat. Phosphorus, instead, must be supplied at concentrations ranging from 37 kg P${}_2$O${}_5$ ha${}^{-1}$ to 77 kg P${}_2$O${}_5$ ha${}^{-1}$ for sunflower, and with concentrations ranging from 56 to 98 kg P${}_2$O${}_5$ ha${}^{-1}$ for wheat (Figure 7).

4. Discussion

Normalized Difference Vegetation Index (NDVI) was estimated from Sentinel 2 satellites constellation data on a 5 ha wheat field crop during the phenological time-step of the start of grain filling period to test whether it can be a reliable proxy for the mineral status of the soil. NDVI from Sentinel 2 is commonly used for crop yield mapping, fertilizer use, and minimizing nitrogen loss to water [8], usually in precision agriculture settings [59] despite its medium-to-low spatial resolution. A limitation of this study concerns the saturation effect that may affect NDVI by losing its linear relationships with aboveground biomass at higher biomass values. The saturation effects may hinder its estimation performance and confound its relationship to soil chemical/physical properties. A multi-sensor approach might help overcoming this effect (see, e.g., in [60]).

Soil mineral status was sampled in the layer where the great part of wheat roots are distributed [37], particularly in clayey, unstructured soils, with high bulk density and worked at shallow depths in the Mediterranean climate, characterized by low rainfall and low soil moisture content [61].

On the one hand, despite the fact that the availability of the three macronutrients (nitrogen, phosphorus, and potassium) is commonly associated to soil fertility and crop growth, they did not account for the spatial variation of NDVI over the study field nor did any of the micronutrients sampled or any other soil physical features. On the other hand, NDVI variability was associated with the soil colloidal status and, in particular, with the components most active in determining the flocculation of the clays and the aggregation of the soil particles, i.e., the organic matter and the Ca${}^{2+}$ ions adsorbed on the exchange complex [62].

The role of SOM is highlighted by the positive correlation between NDVI and the colloids index: even small variations in organic matter can significantly influence the structure of the soil, inducing improvements in soil physical-chemical fertility and in plant nutrition. The presence of higher quantities of Ca${}^{2+}$ ions induces the formation of a larger number of bonds between the mineralogical component and the organic matter, which causes the formation of a higher number of stable soil aggregates [63].

An important ecosystem property that contributes to soil organic carbon (SOC) stabilization and soil structure stability is the interaction between SOC and cations or minerals. SOC can be stabilized by organo-cation or organo-mineral interactions [64]. When the polyvalent cations concentration is high, it becomes sufficient to flocculate and precipitate soluble organic matter. In particular, research in Ca-rich field environments has highlighted a positive correlation between exchangeable Ca${}^{2+}$ and SOC concentration. Ca is a plant macronutrient, and it has a localized positive effect on net primary productivity and soil organic matter inputs both for aboveground and belowground biomass. Exchangeable Ca concentration is correlated with a reduction of SOC leaching, photo-oxidation and respiration [65]. Furthermore, the role of clays can also have different effects on the stabilization of soil aggregates, depending on clay mineralogy, particularly at large clay contents. Wuddivira and Camps-Roach [66], by treating a clayey-kaolinitic soil and a sandy-kaolinitic soil with Ca${}^{2+}$ and organic matter, improved aggregation within a short time, while the same treatment on a clayey-smectic soil gave the opposite effect, suggesting the need for adequate time for aggregate improvement through Ca${}^{2+}$ bridging.

The lack of relation between N, P, and K content in the soil and NDVI has been validated by elaborating fertilization plans both for a successive sunflower crop and for a successive wheat crop elaborated following the current regulations enacted by the competent regional administration. An excess of nitrogen and potassium (${\mathrm{K}}_2\mathrm{O}$) and a slight demand for phosphorus (${\mathrm{P}}_2{\mathrm{O}}_5$) were highlighted by both fertilization plans, although with spatial variations within the field. However, it must be stressed that most of the phosphorus added to the soil is likely to be immobilized, due to the high active limestone content of the soil.

Often, spatial variability of NDVI is commonly associated with nitrogen demand by the crop so that N fertilization plans are deployed by thresholding NDVI into spatially explicit classes and assigning them different N concentrations [17,67]. Although limited to the specific condition of the study, a lower demand for macronutrients may be a less frequent condition than one might think and hence to be worthy of investigation. Should our results be confirmed on wider soil contexts, in Mediterranean intensively used soils the traditional macronutrient fertilization practices could be limited. As the vegetation model has suggested that NDVI variability is to be associated with the variability in SOM, an organic fertilization could be better suited to increase soil matter and crop yield than classic N or NP, or NPK fertilization. Higher SOM tends to mean a larger soil microbial population and therefore potentially higher N supply through mineralization [68]. Second, it should be emphasized that careful evaluation of soil chemical and physical properties should be instrumental to the deployment of properly conceived fertilizations treatments even in precision agriculture frameworks.

5. Conclusions

The commonly used Normalized Difference Vegetation Index from the Sentinel-2 satellite constellation has demonstrated to be very sensible even to the narrow crop productivity variations in the field. For the specific conditions of the study, NDVI variability was influenced by the colloidal status of the soil more than its nutrient availability. Further research is needed to confirm whether the relationship between NDVI and colloid index reported still holds in other clayey soils and in other soil contexts.

Fertilization plans that do not take into consideration soil chemical and physical features may wrongly supply one or more macronutrients under the simplifying assumption that NDVI is solely correlated to nitrogen availability. Variations in crop productivity can be associated with different functional qualities of the soil. Nitrogen, due to its dynamism and its mobility, is the main factor of crop production variability in soils with low fertility. In clay soils, functional qualities can be connected to other factors, first of all the content of organic matter, exchangeable cations, and the quantity and composition of the clay minerals.

Soil and fertility degradation caused by the extensive use of pesticides and fertilizers can be tackled by applying a reasoned analysis of soil properties viewed as a whole system. An integrated agro-ecological assessment coupled to remote sensing can provide useful insights into a more sustainable and targeted fertilization approach. This is particularly true in intensively used soils, such as those the study was based on, where soil organic matter is steadily declining. This approach (i.e., use the NDVI/soil characteristics association) can play a remarkable role to better target the nutrient inputs and to avoid unjustified use of fertilizer.