Assessing the Impact of Water Salinization Stress on Biomass Yield of Cardoon Bio-Energetic Crops through Remote Sensing Techniques

Abstract

Various species of cultivated thistle, such as Cynara cardunculus L. (cardoon), exhibit interesting features for industrial biomass production as bioenergy crops, given also their advantageous adaptation capacities to typical Mediterranean climate trends, with noticeable resilience to drought and salinization stresses. The in situ hyperspectral reflectance responses of three genotypes of cardoon plants, irrigated with water at different salinity levels, have been tested for assessing the effects on their biophysical parameters, aiming at improving the biomass yield for bioenergy production, minimizing at same time the environmental impacts and the exploitation of soils and waters resources. The leaf and canopy reflectance hyperspectral signatures, acquired at three different growth stages with biometric measurements, were statistically analyzed (ANOVA, Tukey’s test, graphs), as noise-resilient spectral indices, sensible to different plant features of interest. Their broadband versions, based on the Landsat 8 OLI and Sentinel 2 MSI satellite sensors, were also evaluated in perspective of operative and extensive remote crop monitoring from space. The results highlighted significant differences in some spectral index responses, related to different cardoon genotypes and water salt concentration. The biometric data supported by red-edge indices modelling evidenced the impact of the highest salt water concentration (200 mM/L) on the plant growth and yield.

1. Introduction

In the perspective of exploitation increase in renewable energies [1,2], various productive crops are currently considered as bioenergy resources capable of supporting reduction in the world’s dependence on fossil fuels and decrease in gas emissions, as already required in Kyoto agreement [3,4,5]. Their relevant features are high biomass production, high calorific value, low agronomic requirements, low cost of production, and fertilizer needs [6]. The currently available advanced precision farming techniques are able to adequately support intensive bioenergy crops production for yield and profit maximization, lowering at same time both water, fertilizer, and environment impact using extensive nondestructive crop monitoring methods. Biomass produced by bioenergy crops can be used for the production of heat or electricity by means of direct combustion or production of biofuel and biogas using pyrolysis and gasification technologies [5]. Some bioenergy crops are suitable for cultivation, especially in marginal lands affected by industrial contamination or salinization, without irrigation to avoid rise in production costs. In this context, cardoon has been found as a low-input annual crop tolerant to drought and various environmental and climate stresses and well-adapted to be exploited as bioenergy crop in the Mediterranean areas [7,8]. Exploitation of the most performing cardoon genotypes in terms of stress tolerance, water deficit, and salinization resilience, with maximization of biomass production and other valuable production outcomes, was the goal of various research [9].

The noninvasive and nondestructive passive proximal sensing technologies can provide effective field tools for quick and early diagnostics in supporting the combat the biotic and abiotic stress factors affecting the plants health as well as the crop production and quality [10,11,12,13]. The optimization of these aspects is the goal of the modern precision farming to face the growing needs in terms of yield, avoiding or mitigating at same time the losses, the damages, and the environmental impact [14,15,16]. In particular, the diagnostic techniques based on the proximal hyperspectral reflectance spectroscopy are spreading significantly in key research fields of applied plant science. These nondestructive methodologies have been widely used also as ground calibration in the growing applications of the crop extensive monitoring, using remote sensing approaches based on aerial and satellite multi/hyperspectral sensors [17,18,19]. In particular, due to continuous and repetitive monitoring capability coupled with recent advances in sensors and platform, the satellite remote sensing applications currently provide a broad spectrum of effective and tailored solutions for terrestrial and marine vegetation assessment also linked to related environmental impacts and water quality parameters [20,21,22,23,24]. In this context, the present study has aimed at investigating whether the spectral reflectance responses of three different genotypes of cardoon, subjected to different salinity levels during their first growth stages of Biologische Bundesanstalt, Bundessortenamt, CHemischeIndustrie (BBCH) scale [25,26], could be used to suitably estimate the effects on different plant biophysical parameters (i.e., Leaf Area Index (LAI) and leaf biomass density). The reflectance spectroscopy exploiting the solar spectrum in the Uv-Vis-NIR-SWIR ranges was tested as a rapid screening method for monitoring differentiation of cardoon genotypes responses to salt stress. According to the current radiative modelling approaches of vegetation, where the leaf models are nested within those devoted to canopy simulation, different hyperspectral signatures at leaf and canopy levels were acquired. The goal was to capture responses more directly linked, respectively, with plant tissue leaf features/optically active substances concentration and plant 3D structure/LAI and related biomass distribution [27]. The collected spectral signatures were then preprocessed and exploited, singularly and under the form of noise-resilient vegetation spectral indices sensible to different plant features and in agreement with those implemented in literature and previous works [28,29,30]. These spectral indices were profitably used for the retrieval of vegetation biophysical parameters from its hyperspectral signatures, taken from a close distance from the canopy minimizing the noises arising from soil, light diffusion, atmospheric conditions, and the sensor observation [31,32]. However, in general, the observed areas and crop portion are very limited, going from single plant to small group of the growing plants in few square meters/centimeters of ground area observed for each measure. Thus, to provide more effective monitoring capability of these bioenergetic crops, the last generation of high resolution (HR) sensors on board of Landsat 8 OLI (Operative Line Images) and Sentinel 2 MSI (MultiSpectral Instrument) satellites have been evaluated considering their suitable ground resolution (30 and 10 m, respectively) and enhanced spectral features [14,33]. Thanks to their systematic acquisitions of visible (VIS), NIR, and SWIR reflectance data over the entire Earth’s surface, they can provide a more effective support to operational monitoring of growing crops and agri-ecosystems, taking into account the management needs of precision agriculture including yield optimization, and minimizing at the same time the environmental impact, according to cultivation seasonal cycles [16,29]. From this perspective, the multispectral data collected by the Landsat 8 OLI and Sentinel 2 satellite sensors were also considered for extensive and operational monitoring of the salinization effects on cardoon growth. To this end, the in situ acquired hyperspectral signatures were transformed into broadband responses, compatible with the spectral configuration of these satellite sensors, in order to simulate the reflectance data acquired by them. Then, the related spectral broadband indices were analyzed in order to evaluate their capability to capture the salinization effects on the various cardoon genotypes tested here, even from space. From the results, a significant difference of some spectral indices derived from pre-processed narrow and broadband responses of different cardoon genotypes was evidenced by ANOVA analysis, performed on the whole dataset including the three campaigns data. Various spectral indices, such as the widely used Normalized Difference Vegetation Index (NDVI) and the red-edge-based ones, demonstrated a different response to water salinity levels with an additional capability of some of them to discriminate the diverse responses to water salinization of the cardoon genotypes used in the experiment. The description of expressions, acronyms, and abbreviations used in this paper is reported in

Table A1. Abbreviation, expressions and acronyms.

Expression	Description
(Adjusted) R square	(Adjusted) R2 correlation coefficient
a.g.r.	Above ground resolution
ANOVA	Analysis of Variance
BBCH	Biologische Bundesanstalt, Bundessortenamt, CHemischeIndustrie
Bio	Leaf biomass density (mg/cm2)
bx	Broadband spectral response in band/spectral range x
cmpn	Measurement campaign (n:1, 2, 3)
DF	Degree of Freedom, i.e. number of measurements minus parameters
Fact.	ANOVA factor
F Value	ANOVA F statistic
fij	Broad band j filter function
gen	Cardoon genotype
LAI	Leaf Area index
mM	10−3 moles
MSI	MultiSpectral Instrument
NIR	Near Infrared
OLI	Operational Land Imager
plat	plateau
Pl X	Plateau at different water salinization levels (X: A, B, C)
Prob > F	p-value for the whole regression model
p-value	Significance level of single factor/global relationship model
R	R Spearman correlation coefficient
redg	Red Edge
rx	Narrow band spectral reflectance in band x
SNR	Signal to noise ratio
SWIR	Short Wave Infrared
TIR	Thermal Infrared
Typ.	Typology
Uv	Ultraviolet
VIS	Visible
Y	Biomass yield (t/ha)

of the Appendix A.

2. Materials and Methods

2.1. Experimental Setup

Three domesticated cardoon genotypes, whose seeds have been provided by ENEA germplasm collection, were grown in hydroponics boxes located in a greenhouse of the ENEA’s Casaccia research center [25], until the first growth stages of BBCH scale [26], shown in Figure 1.

Three clones of C. cardunculus were previously selected for the morphological characterization on three plants/genotype and subsequently plant selections were performed [26]. The main descriptors applied were plant height and width, length of the youngest well-expanded leaf, stem diameter, number of heads produced, and plant dry weight. The dry matter content of each genotype was determined oven-drying the plant material (heads and biomass) at 103 °C until a constant weight was reached. During their vegetative development in nutrient solution supplemented with different concentrations of NaCl, the plants of the three genotypes were monitored at three different growth stages through hyperspectral reflectance signatures acquisition and calibration by means of leaf biometry parameters. The three cardoon genotypes were grown by means of hydroponic technique in an unheated greenhouse to adequately manage the different salt concentration in their irrigation water. The seeds of the three cardoon genotypes have grown on an inert substrate of standard polystyrene plateau (six rows of 10 holes each (diameter = 4.2 cm) with dimensions of 53.2 × 33.3 × 6.0 cm). The seeds of each specific genotype were planted in the 20 holes of two contiguous rows of the six available, maintaining the same sequence for all the three plateau. Each plateau was then placed in hard plastic boxes of the volume of about 0.05 m³ with a 50 liter water nutritive solution, added with 0, 100, and 200 mM of NaCl, respectively. Crop management of the experimental fields was conducted under low energy inputs (manual weeding, basic fertilization using N 100 kg·ha⁻¹, P₂O₅ 9 kg·ha⁻¹, K₂O 12.5 kg·ha⁻¹). The hydroponic technique of floating system was applied under controlled conditions in the greenhouse. At weekly intervals for the whole duration of the test (35 days), the indoor temperature (18–24 °C), the pH (6.1 to 7.5), and the conductivity (mS·cm⁻²) of the saline solution inside the tanks were monitored. Figure 2 shows the cardoon plantlets growing in the three plateaus at the time of the second spectral acquisition. To keep constant the quantity of the nutritive solution present in the containers, taking into account the absorption of the plants and the evaporation, water with the same NaCl concentration was continuously added, controlling at same time the pH and electrical conductivity. The boxes were placed in the greenhouse with a disposition capable of providing mostly homogeneous lighting of the sun spectrum on all three plateaus. The seeds in the three plateaus were planted at the beginning of March 2016. The first acquisition of spectral signatures was accomplished on 7 April, the second took place on 21 April with contemporary measurements of biomass, Leaf Area Index (LAI), and leaf biomass density (Bio), while the last spectral acquisition was conducted on May 11. Between acquisitions, there was a plant development interval of about 12–14 days. After the third acquisition, the experiment was interrupted when the cardoon plants in the plateaus irrigated with unsalted water reached the 1.5–1.7 stage of BBCH growing scale, corresponding to about 100% of ground cover, while those irrigated with most salinized water were at 1.2–1.4 BBCH stage.

Considering the not excessively large population and the defined geometric configuration of the experiment, based on the characteristics of the artificial substrates used, a systematic sampling scheme was adopted for the acquisition of spectral signatures. This allowed us to exploit the repetitiveness of the measurements on the plants of each genotype present in the three boxes with different salinity of the water during the measurement campaigns. During the three campaigns, the field measurements have been carried out trying to operate with the most favorable and constant (i.e., near noon) sun illumination.

2.2. Hyperspectral Data Acquisition

The hyperspectral reflectance data have been acquired through Advanced Spectral Device (ASD) FieldSpecPro hand-held spectrometer/radiometer that allows one to collect point radiance/reflectance responses between 350 and 2500 nm wavelengths with a spectral resolution ranges of 3 nm (within 350–1000 nm) and 10 nm (within 1000–2500 nm). The measurement probe was equipped with standard foreoptics, which provided field of view by 25°, employed in all acquisition campaigns. To acquire reflectance data modulated by the sunlight interaction with cardoon plants, a normalization with the solar irradiance spectrum collected by a spectralon (polytetrafluoroethylene) surface (characterized by flat reflectance Lambertian response) was used. In the used device, this normalization spectrum is registered once at the beginning of the measurements, assuming constant illumination during the subsequent reflectance data acquisition. In case of changes in solar lighting, acquisition of the normalization spectrum should be repeated in order to avoid errors in reflectance measurements. As absorption effects of water vapor concentration are extremely variable in two ranges around 1380 and 1910 nm, here, the measured reflectance values are very unstable and oscillating between extremes much higher than −1 and 1 also arising from an unreliable calibration, so these data were excluded. Where there were values higher than 1 in the valid wavelength ranges caused by an under-calibration, a renormalization at the highest reflectance was performed.

The acquisition of spectral signatures was carried out with a vertical probe (Figure 3) connected to the radiometer through the optical fiber, perpendicular to the reflecting surfaces of the leaves or over the plants. As previously stated, the water salinization affects not only the concentration of nutrients and the tissues health at the leaf level but also the Leaf Area Index (LAI) of the entire plant, directly connected to the canopy architecture and to the biomass production. In accordance with the usual light–vegetation interaction model acting at the leaf and canopy levels, two different empirical approaches were introduced for spectral measurements to take account for these two nested effects. In particular, taking into account the 25° field of view (FOV) instrument, measurements, the signatures of both single leaves and the entire plant (canopy) were acquired by changing the distance between the perpendicular probe and the reflecting surfaces (Figure 2).

2.3. Hyperspectral Vegetation Indices

In general, due to the peculiarity of the typical spectral signature of the plant photosynthetic apparatus, its single wavelength spectral band responses in different spectral ranges are widely exploited under the form of specific combinations called spectral vegetation indices (SVI). SVI improves the capability to capture the spectral reflectance behaviors linked to different biophysical parameters of the plants, minimizing at same time the noise contributions from various noise factors (i.e., soil reflectance, spectral response anisotropy, atmospheric turbidity). The SVIs derived particularly from red and near infrared (NIR) reflectance signals are widely used for monitoring the development of specific vegetation features and estimating the biophysical and yield traits of agricultural crops (e.g., biomass, water content, pigments). In particular, Normalized Difference Vegetation Index (NDVI) was usefully exploited for assessing various important biophysical parameters of the plants characterizing their productivity and health. These SVI-based methods have been exploited to study forestry, grasslands, and rangelands since the early years of satellite EO [34]. Due to their formulation, the ratio-based SVIs, such as NDVI, have an intrinsic capability to minimize the multiplicative noise associated to variable topography and atmosphere [29], while, at least in their basic form, are less effective in decreasing the influence of variable soil reflectance. For instance, to minimize this latter factor, mainly arising during the various development stages of plants that involve rising percentages of soil covered by green leaf surfaces, specific SVI have been developed like Soil Adjusted Vegetation Index (SAVI) and Optimized Soil Adjusted Vegetation Index (OSAVI), which demonstrated various levels of effectiveness in different situations. Some indices, more sensible to chlorophyll and leaf content of optically active substances, were introduced to better exploit the leaf hyperspectral data. Taking into account the strong link of the red-edge wavelength reflectance with LAI and stress status of vegetation arising from various factors, different indices derived from responses in this spectral range were assessed [27], even under the form of Red-Edge Inflection Point (REIP) of the hyperspectral signature trends [35]. The water indices used are characterized by their response to leaf/plant water content indirectly linked with biomass and LAI. In

Table 1. Hyperspectral vegetation indices. (Typ. = index typology).

Typ.	Index (Denomination)	Formulation	Reference
1	NDVI (Normalized Difference Vegetation Index)	(rNIR − rred)/(rNIR − rred)	[36]
1	SAVI (Soil Adjusted Vegetation Index)	(rNIR − rred) (1 − L)/(rNIR − rred + L)	[37]
1	OSAVI (Optimized Soil Adjusted Vegetation Index)	(rNIR − rred) 1.16/(rNIR − rred + 0.16)	[38]
1	NDVI1 (Normalized Difference Vegetation Index Red Edge)	(r807 − r736)/(r807 + r736)	[13]
2	TCARI (Transformed Chlorophyll Absorption Ratio Index)	3[(r700 − r670) − 0.2(r700 − r550)(r700/r670)]	[39]
2	TVI (Triangular Vegetation Index)	0.5[120(r750 − r550) − 200(r670 − r550)]	[40]
2	SIPI (Structural Insensitive Pigment Index)	(r800 − r445)/(r800 + r680)	[41]
2	NPCI (Normalized Pigment Chlorophyll Index)	(r680 − r430)/(r680 + r430)	[42]
2	MCARI (Modified Chlorophyll Absorption Ratio Index)	1.2[2.5(r800 − r670) − 1.3(r800 − r550)]	[43]
3	RE11 (Red Edge 1)	(r750 − r700)/(r750 − r700)	[44]
3	RE21 (Red Edge 2)	(r740 − r720)/(r740 − r720)	[45]
3	RIRE (Reflectance Index Red Edge)	(r750 − r800)/(r69 5+ r740)	[46]
3	REIP1 (Red Edge Inflection point)	700 + 40{[(r670 + r780)/2] − r700}/(r740 − r700)	[44]
3	REIP2 (Red-Edge Inflection point)	705 + 35{[(r665 + r783)/2] − r705}/(r740 − r705)	[44]
4	NDWI (Normalized Difference Water Index)	r860 − r1240/r860 + r1240	[47]
4	NDWI1 (Normalized Difference Water Index)	r858 − r1640/r858 + r1640	[48]

, the exploited indices derived from the hyperspectral signatures collected during the three monitoring campaigns were reported; the generic reflectance response is indicated with r and the reference wavelength as subscript.

The first three of selected indices are based on the spectral characteristic of vegetation absorbing the red wavelengths for its photosynthetic primary function and with a reflectance maximum in the NIR ranges due to optical properties of its vegetal tissue strata. They are more sensible to plant structural features and canopy biophysical parameters, such as LAI and biomass, while the other ones are formulated to capture mainly the responses of optically active substances at leaf level, such as chlorophyll or pigments. In these narrow band indices, the NIR and red wavelengths were assumed, respectively, at 800 and 670 nm. The other SVIs are related to so-called red-edge spectral range, located between red and NIR wavelengths, where the spectral reflectance of photosynthetic vegetation exhibits a sharp rise and concavity change with inflection point. These red-edge-based indices have proven effective in assessing LAI and detecting the stress status, with enhanced robustness against the structural variability of plants [49]. The spectral responses of last two spectral indices are related mainly to water absorption bands in the SWIR ranges and make them sensible to water leaf content. According to specific interactions with vegetation, the selected indices exploit spectral responses extending across a broad wavelength intervals going from visible to NIR and SWIR of the acquired hyperspectral signatures. In synthesis, according to their specific features, the selected indices can be included in the four typologies: structural/LAI (Typ. 1); chlorophyll/pigments (Typ. 2); red edge (Typ. 3); water content (Typ. 4).

2.4. Biometry Data

In the second campaign, after the hyperspectral signatures, from the cardoon plants with a sufficient development degree, nine samples for each plateau were collected to exploit them for laboratory estimate of LAI and biomass biophysical parameters. To minimize the disturbance to the canopy spectral response, the nine plantlets were taken from the edge of each plateau, three for each genotype, for a total of 27 samples. Subsequently, the leaves of each plant were counted and their length and maximum width were measured. Then, their weight (g) was also derived to assess the means for each sample. The average leaf surface for each plant was calculated assuming a fixed ratio of 0.65 from the rectangular shape, then the total surface was obtained multiplying it for the number of leaves. The assessment of the plant LAI was accomplished using the circular surface corresponding to about 12 cm of diameter of the radiometer FOV for canopy spectral measurement as normalization factor. This diameter is equivalent to the width of the two row holes of the plateau used for the plant of the same cardoon genotype. Finally, the biomass and LAI averages for each genotype and plateau were assessed in order to provide a reference for the other spectral data.

The cardoon yield in terms of aboveground biomass for hectare Y (q/ha) is one of the most important parameters for bio-refinery and bioenergy production, whose assessment and maximization is the primary interest of this research. The three measurement campaigns were carried out in the first growth stages of the cardoon plantlets, when they have leaves only without stalk and heads that, in the field, develop subsequently during the last development periods, before drying and harvesting. Thus, in this lab condition, the fresh aboveground biomass potentially provided, approximately can be estimated from Bio and LAI measured parameters, to assess the productivity Y without considering the stalks and heads contribution. Although, in general, Y has essentially an agronomic meaning and is not well sounding in this experimental controlled environment where cardoon plants were grown using hydroponic technique without any fertilizer, it was, however, introduced as follows in agreement with the perspective of extensive field assessing of the bioenergy crops in the framework of precision farming approaches:

Y (q/ha) = Bio (mg/cm²) × AI (m²/m²)

(1)

2.5. Broadband Indices

The Landsat 8 OLI (L8-OLI) and Sentinel 2 MSI (S2) data consist of systematic multispectral measurements of the Earth’s surface reflectance in broad spectral bands (ranging between visible (VIS), near infrared (NIR), and short-wave infrared (SWIR)/thermal infrared (TIR) wavelength), taken by the satellite sensors, orbiting at 700–800 km far from surface. The close-range operation of the ground spectroscopy allows the acquisition of spectral signatures sampled at many narrow bands (2–10 nm wide), while to make the remote sensor capable of capturing sufficient reflected energy from Earth’s surface, few broad spectral bands (20–200 nm wide) acquisition are usually exploited in the remote sensing techniques. These high-resolution (HR) sensors feature more bands, an improved signal-to-noise ratio (SNR) with a more effective radiometry (12 bit), enhanced spatial resolution (10–30 m), and revisit capabilities, better than the older Landsat HR sensors, which have been exploited for early application in agriculture since the start of space Earth-observation (EO) activities. At present, the multispectral data systematically provided by the dual satellite systems S2 with its 10 m of above ground resolution (a.g.r.), constitute the basic information for ESA (European Space Agency)-E.U. Copernicus operational services, in many strategic sectors, in particular to preserve the environment and sustain productivity in agriculture by assessing agricultural land use and trends, crop health conditions and yield forecasts. The applications of Copernicus services, in addition to precision farming, include seasonal mapping of cultivated areas, water management, and drought monitoring, as well as subsidy controls.

The L8-OLI (11 bands) and S2 (12 bands) sensors use electronic bandpass filters to derive the spectral responses for the various acquisition broadbands from the full spectrum radiance (Figure 4), captured by their optical telescopes (

Table A2. Landsat 8 OLI acquisition bands. The band width is evaluated as FWHM (Full Width Half Maximum) of the related filter function.

Bands	Landsat 8 OLI Band	Center Wavelength (nm)	Band Width (nm)	Spatial Resolution (a.g.r.) m
b1	Coastal	442.96	15.98	30
b2	Blue	482.04	60.04	30
b3	Green	561.41	57.33	30
b4	Red	654.59	37.47	30
b5	NIR	864.67	28.25	30
b6	SWIR1	1608.86	84.72	30
b7	SWIR2	2200.73	186.66	30
b8	Pan	589.50	172.40	15
b9	Cirrus	1373.43	20.39	30

and

Table A3. Sentinel 2 MSI acquisition bands. The band width is evaluated as FWHM (Full Width Half Maximum) of the related filter function.

Bands	Sentinel-2 MSI Bands	Central Wavelength (nm)	Band Width (nm)	Spatial Resolution (a.g.r. m)
b1	Band 1: Coastal	443	20	60
b2	Band 2: Blue	490	65	10
b3	Band 3: Green	560	35	10
b4	Band 4: Red	665	30	10
b5	Band 5: Red-edge 1	705	15	20
b6	Band 6: Red-edge 2	740	15	20
b7	Band 7: Red-edge 3	783	20	20
b8	Band 8: NIR	842	115	10
b9	Band 8A: NIR narrow	865	20	20
b10	Band 9: Water Vapor	945	20	60
b11	Band 10:SWIR Cirrus	1375	30	60
b12	Band 11: SWIR	1610	90	20
b13	Band 12: SWIR	2190	180	20

of Appendix A). These filters are characterized by analytical filter functions (f_ij) made available by the data provider (the graphs showing the spectral bands filters functions of L8-OLI and S2 MSI sensors are reported in Figure A1 of complementary materials). These filter functions have been exploited to obtain the spectral responses compatible with those of satellite broadbands, from the hyperspectral signatures (r_i) previously acquired on field, by means of convolution:

$${\mathrm{b}}_{\mathrm{j}}={\displaystyle \sum }_{\mathrm{i}}{\mathrm{f}}_{\mathrm{i},\mathrm{j}}{\mathrm{r}}_{\mathrm{i}}/{{\displaystyle \sum }}_{\mathrm{i}}{\mathrm{f}}_{\mathrm{i},\mathrm{j}}$$

(2)

where b_j states for broadband response in band j, f_ij is the filter function of j broadband of satellite sensor and r_i the hyperspectral reflectance data referring to i narrow band. To obtain the various spectral band responses, it is required to use the different f_ij of the specific sensor, according to the previous approach adopted for hyperspectral data; also, the broadband responses have been exploited under the form of specific vegetation spectral indices widely, used for assessment of various biophysical parameters of vegetation and crops [50,51]. The following

Table 2. Landsat 8 OLI and Sentinel 2 MSI broadband indices.

Typ.	Index	Formulation	Reference
1	NDVIL	(bNIR − bRed)/(bNIR + bRed)	[34]
1	SAVIL	(1 − L)(bNIR − bRed)/(bNIR + bRed + L)	[35]
1	OSAVIL	(1+ 0.16) (bNIR − bRed)/(bNIR + bRed + 0.16)	[36]
1	EVIL	2.5 (bNIR − bRed)/(bNIR + 6 bRed − 7.5 rBlue + 1)	[51]
2	SIPIL	(bNIR − bCoastal)/(bNIR + bRed)	[39]
2	NPCIL	(bRed − bCoastal)/(bRed + bCoastal)	[40]
4	NWDI1L	(bNIR − bSWIR1)/(bNIR − bSWIR1)	[45]
4	NWDI2L	(bNIR − bSWIR2)/(bNIR − bSWIR2)	[45]
3	RE11s	(bredg2 − bredg1)/(bredg2 + bredg1)	[53]
3	RE21s	(bredg3 − bredg1)/(bredg3 − bredg1)	[53]
3	REIPs	740 + 40 × (bred + bredg3)/(2 − bredg1)/(bredg2 − bredg1)	[54]

includes the spectral vegetation indices assessed from the field hyperspectral signatures using the L8 OLI filter functions (subscript L). Analog indices have been estimated using S2 spectral band filters (subscript s) functions of sensor focusing on its red edge and different NIR acquisition channels, some of them were reported in Table 2 [52]. In particular, for S2 different formulations of NDVI have been introduced, using, respectively, its NIR, red edge, and water vapor bands. Three red-edge indices have been tested, the first two (RE11s, RE21s) are under form of normalized difference index of the three red-edge acquisition bands [53], exploited also as Red-Edge Inflection Point (REIP). Finally, additional indices sensible to vegetation water content similar to those of OLI sensor have been assessed by means of narrow NIR and two SWIR S2 responses.

2.6. Statistical Analysis and Models

The ANOVA analysis of spectral indices derived from both the leaf and plant canopy hyperspectral signatures acquired during the three field campaigns has been performed. Two-way ANOVA was carried out assuming both the water salinization level of the three growing plateaus and the cardoon genotypes as two variation factors. The same approach was used for the Landsat 8 OLI and Sentinel 2 MSI broadband indices assessed as described above from the hyperspectral data acquired at canopy level. In this way, all the spectral indices including the broadbands ones were tested for their discrimination capability of differences between the reflectance of the cardoon plants subjected to different stress from the three salinization levels of irrigation water. In addition, an assessment of their further aptitude to distinguish between the three cardoon genotypes was provided. From the ANOVA results, the most performing indices at leaf and canopy levels were further evidenced through Tukey’s test to refine the characterization of their detection capability based on comparing the results for the pairs derived from the combination of the different water salinization levels and genotypes. This statistic test allowed to evaluate the global performance of each spectral index utilized. In particular, the test assesses the number of pair combinations discriminated for spectral effects of the three different water salinization levels (plat columns) and genotypes (gen columns) with p-value < 0.05 (n > 0.05). In the cells of the results tables, the sum (with a maximum of three) of the pairs discriminated was reported as a proxy of global performance of the index, respectively, for the water’s salinization level and genotype. In case of null sum, the cell is left empty.

The biometry data, as LAI (m²/m²) and Bio (mg/cm²) parameters assessed from the cardoon plants samples collected during the second campaign of the 21-04-2016 with the related hyperspectral signatures, were firstly analyzed through the two-way ANOVA approach to evidence their dependence from the water salinization levels and cardoon genotypes. In order to allow the modelling through the above introduced spectral indices averages of these biophysical parameters, estimates for the three different water salinization levels and genotypes were exploited. Considering that the LAI and the Y (Equation (1)) are more directly linked with biomass resources production by bioenergy crops, their modelling capability through spectral canopy indices, including the broadband ones, was estimated as Spearman’s correlation coefficient. The latter ranges between −1 and 1 and allowed us to identify the most performing indices through an effective and comparative graphical approach as a bidimensional scatter plots. Subsequently, the linear models based on the selected indices were assessed and analyzed.

3. Results and Discussion

3.1. Field Data

Table 3. Number of different field measurements for each campaign (cmp1; cmp2; cmp3).

Typ.	cmp1	cmp2	cmp3	Total
	7 April 2016	21 April 2016	11 May 2016
Leaf	90	36	27	153
Canopy	27	9	27	207
Biometry		27		27
total	117	72	54	243

includes the number and type of field measurements acquired for each campaign on the cardoon plantlets of the three plateaus. In particular, in the cmp1 campaign, 90 leaf spectral signatures were collected, 30 for each plateau containing water at different salt concentration, and divided equally for the three genotypes. In the same campaign, nine canopy spectral measurements for each plateau, three for each genotype, were carried out. In the second campaign, 12 leaf spectral signatures (four for each genotype) and three canopy spectral data (one for each genotype) for each plateau were collected. These field data included biometrics measurements, performed in laboratory for each plateau and genotype using samples collected during spectral acquisitions. In the last campaign (cmp3), three spectral signatures for the plants of each genotype growing in the three plateaus at different water salt concentration were finally collected. A total of 360 leaf and canopy spectral signatures (of three campaigns) with laboratory assessment of 27 biometric values (cmp2 campaign) of biophysical parameters (specific biomass and leaves surface) related to the three genotypes for each salinization level were provided for subsequent statistical analysis.

3.2. Hyperspectral and Broadband Indices

The graphs of Figure 5 show the canopy hyperspectral signatures averaged for each of the three cardoon genotypes growing in the plateaus A, B, C, with rising water salinization concentrations, for all three campaigns. The similar graphs obtained from leaf spectral signatures were reported in the auxiliary materials (Figure A2). The graphs along the column show the hyperspectral signatures of cardoon genotypes growing at different salinization degrees in the three campaigns. All graphs highlight the instabilities of the reflectance values recorded around the 1380 and 1890 nm wavelength ranges corresponding to air water vapor effects, with reflectance rapid variations exceeding the physical threshold (0, 1). In addition, as shown in the graph of the first campaign (cmp1), for plateau A (pl A), the effects due to the calibration fault led to unreliable reflectance values higher than 1. The spectral signatures of the different genotypes are sufficiently distinguishable in all the graphs, while those referred to plateaus A and B demonstrated a general decrease in spectral VIS responses and increase in NIR ones in the subsequent campaign, according to the increase in the chlorophyll absorption and leaves number of developing plantlets. The plants developed with the highest water salinization (pl C) showed instead an inverse distribution of the genotype hyperspectral signatures in the first (cmp1) and in the last (cmp3) measurement with an agreement on their trends of intermediate campaign (cmp2). In the first campaign, the development of most plants is in the early stages, and due to their low LAI values, the background reflectance contribution was significant. In the second campaign, the hyperspectral reflectance curves of the different plant genotypes stabilized and increased their differentiation in the third one.

Although most of the undesirable effects of noise above evidenced were corrected and mitigated in a pre-processing step, the exploitation of hyperspectral signatures under the form of normalized ratio spectral indices allowed us to further reduce this problem.

All the leaf and canopy (narrow band) indices introduced above were derived from the pre-processed hyperspectral data, then also the broadbands ones were assessed using the Landsat 8 OLI and Sentinel 2 MSI spectral bands’ filter functions (Figure A1). The results obtained from the two-factor ANOVA analysis, through Tukey’s test, were summarized in

Table 4. ANOVA Tukey’s test analysis results summary obtained for leaf and canopy hyperspectral indices related to the three measurement campaigns (cmp1, cmp2, cmp3). The plat and gen labels indicate, respectively, the plateau (different water salinization) and cardoon genotype factors.

Indices	Leaf Reflectance						Canopy Reflectance						Total
	cmp 1		cmp2		cmp3		cmp 1		cmp2		cmp3
	Plat	Gen	Plat	Gen	Plat	Gen	Plat	Gen	Plat	Gen	Plat	Gen
NDVI	2	1	2		1		2				1		9
SAVI	2		1				3		1		2	2	11
OSAVI	1		2				2		1		3	1	10
NDVI1	3	1	2				2	1			1		10
TCARI	1		3				2		1		3	2	12
TVI		1	2								2	2	7
SIPI	2	2	2		1		2				1		10
NPCI	2		2		1		2				2	1	10
MCARI			2				2		1		1	1	7
RE11	3		2				3				1		9
RE21	2		2		1		3		1		1		10
RIRE	2	1			3		2				1		9
REIP1	2		2				2		1	1	1		9
REIP2	2		2				2		2	1	1		10
NDWI	3	2	2		2		2			1	1		13
NDWI1	2	2			1		2			1	2	1	11

. For each of the selected spectral indices, the table includes the number of combination pairs referred to the plateau (plat. label) and genotype (gen. label) found significantly different (discriminated) in terms of means (p-values < 0.05).

Table 5. ANOVA Tukey’s test analysis results summary obtained for Landsat 8 OLI broadband canopy indices (the subscript _L states for Landsat 8 OLI sensor) related to the three campaigns (cmp1, cmp2, cmp3). The plat and gen labels indicate, respectively, the plateau (different water salinization) and cardoon genotype factors.

Indices	cmp 1		cmp2		cmp3		Total
	Plat	Gen	Plat	Gen	Plat	Gen
NDVIL	2				2		4
SAVIL	2		2		2		6
OSAVIL	2		1		2		5
EVIL	3		2		1		6
SIPIL	2				2		4
NPCIL	2				3	1	6
NWDIL	2				2	1	5
NWI1L	2			1	2	1	6

and

Table 6. ANOVA Tukey’s test analysis results summary obtained for S2 MSI broadband canopy indices (the subscript s states for S2 sensor) related to the three campaigns (cmp1, cmp2, cmp3). The plat and gen labels indicate, respectively, the plateau (different water salinization) and cardoon genotype factors.

Indices	cmp 1		cmp2		cmp3		Total
	Plat	Gen	Plat	Gen	Plat	Gen
NDVI1s	3				2		5
NDVI2s	2				3		5
NDVI3s	3				2		5
SAVIs	2		2		2		6
OSAVIs	2		1		2		5
EVIs	2		2		1		5
SIPIs	2				3		5
NPCIs	3				2	1	6
RE11s	2				2		4
RE21s					2		2
NDVI4s	3				2		5
REIP1S	2		1	1	2		6
REIP2S	2		1	1	2		6
REIP3S	2		1	1	2		6
NWDIs	2				3		5
NDWI1s	2			1	2	1	6
NDWI2s	2				2		4

show instead the number of pairs discriminated by the broadband indices (at canopy level) assessed, respectively, for the L8 OLI and S2 MSI satellite sensors. The corresponding ANOVA detailed p-value tables for all the indices were provided within the Appendix A (

Table A4. p-values obtained from two way ANOVA analyses of structural/LAI and Chlorophyll/pigments leaf hyperspectral indices. The values lower than 0.05 have been reported in bold.

cmp	Fact.	NDVI	SAVI	OSAVI	TCARI	TVI	SIPI	NPCI	MCARI1
1	gen	6.19 × 10−4	7.51 × 10−1	1.43 × 10−1	9.61 × 10−1	5.52 × 10−3	1.28 × 10−4	1.19 × 10−1	7.47 × 10−1
	plat	3.24 × 10−5	2.87 × 10−7	1.85 × 10−4	2.92 × 10−3	9.98 × 10−1	2.31 × 10−3	2.17 × 10−5	3.51 × 10−2
	model	1.33 × 10−6	3.11 × 10−6	3.65 × 10−4	1.83 × 10−2	3.25 × 10−2	1.15 × 10−5	5.86 × 10−5	1.16 × 10−1
2	gen	2.79 × 10−1	3.72 × 10−1	3.59 × 10−1	4.77 × 10−1	3.63 × 10−1	2.63 × 10−1	5.42 × 10−1	3.77 × 10−1
	plat	1.60 × 10−5	2.60 × 10−2	4.80 × 10−3	2.40 × 10−6	1.15 × 10−3	1.23 × 10−5	1.95 × 10−5	3.70 × 10−4
	model	7.77 × 10−5	5.76 × 10−2	1.39 × 10−2	1.72 × 10−5	4.10 × 10−3	5.93 × 10−5	1.28 × 10−4	1.54 × 10−3
3	gen	1.49 × 10−1	1.82 × 10−1	1.35 × 10−1	1.78 × 10−1	1.49 × 10−1	1.30 × 10−1	9.81 × 10−1	1.64 × 10−1
	plat	4.44 × 10−3	2.47 × 10−1	1.52 × 10−1	3.63 × 10−1	1.75 × 10−1	4.52 × 10−3	3.42 × 10−2	9.43 × 10−2
	model	7.98 × 10−3	1.93 × 10−1	1.12 × 10−1	2.47 × 10−1	1.33 × 10−1	7.45 × 10−3	1.32 × 10−1	9.08 × 10−2

,

Table A5. p-values obtained from two way ANOVA analyses of Red-Edge and water content leaf hyperspectral indices. The values lower than 0.05 have been reported in bold.

cmp	Fact.	RE1	RE2	RE11	RE21	RIRE	NDVI1	REIP1	REIP2	NDWI	NDWI2
1	gen	6.08 × 10−1	6.67 × 10−1	2.95 × 10−1	4.67 × 10−1	1.16 × 10−3	1.35 × 10−2	8.02 × 10−1	7.89 × 10−1	4.01 × 10−4	1.49 × 10−5
	plat	4.53 × 10−4	7.12 × 10−8	6.25 × 10−11	2.45 × 10−13	4.80 × 10−11	3.75 × 10−13	7.15 × 10−12	7.23 × 10−12	0.00	1.96 × 10−11
	model	2.35 × 10−3	7.45 × 10−7	6.12 × 10−10	3.94 × 10−12	1.75 × 10−11	9.32 × 10−13	1.33 × 10−10	1.34 × 10−10	0.00	5.42 × 10−13
2	gen	3.78 × 10−1	5.13 × 10−1	2.09 × 10−1	1.26 × 10−1	2.04 × 10−1	9.36 × 10−2	2.17 × 10−1	2.17 × 10−1	4.48 × 10−1	2.29 × 10−1
	plat	1.48 × 10−2	4.94 × 10−2	1.08 × 10−6	4.84 × 10−8	6.74 × 10−1	1.71 × 10−4	1.55 × 10−6	1.55 × 10−6	2.29 × 10−4	1.22 × 10−1
	model	3.68 × 10−2	1.19 × 10−1	5.76 × 10−6	2.68 × 10−7	4.05 × 10−1	3.56 × 104	8.17 × 10−6	8.17 × 10−6	1.10 × 10−3	1.35 × 10−1
3	gen	1.55 × 10−1	1.76 × 10−1	1.24 × 10−1	4.46 × 10−2	3.14 × 10−5	1.61 × 10−1	9.18 × 10−1	9.17 × 10−1	1.38 × 10−11	3.19 × 10−2
	plat	1.40 × 10−1	2.21 × 10−1	7.72 × 10−1	8.89 × 10−1	9.72 × 10−1	8.48 × 10−1	5.98 × 10−1	5.17 × 10−1	1.44 × 10−1	2.13 × 10−1
	model	1.16 × 10−1	1.75 × 10−1	3.08 × 10−1	1.54 × 10−1	2.39 × 10−4	3.92 × 10−1	8.71 × 10−1	8.18 × 10−1	1.21 × 10−10	4.80 × 10−2

,

Table A6. p-values obtained from two way ANOVA analyses of structural/LAI and Chlorophyll/pigments canopy hyperspectral indices. The values lower than 0.05 have been reported in bold.

cmp	Fact.	NDVI	SAVI	OSAVI	TCARI	TVI	SIPI	NPCI	MCARI1
1	plat	3.79 × 10−6	3.94 × 10−7	1.27 × 10−6	6.55 × 10−8	1.00	2.69 × 10−7	6.80 × 10−4	4.62 × 10−8
	gen	1.30 × 10−1	2.33 × 10−1	1.85 × 10−1	4.77 × 10−1	2.41 × 10−1	1.10 × 10−1	2.14 × 10−1	4.22 × 10−1
	model	1.64 × 10−5	2.44 × 10−6	6.78 × 10−6	5.30 × 10−7	5.62 × 10−1	1.38 × 10−6	2.02 × 10−3	3.69 × 10−7
2	plat	1.35 × 10−1	1.23 × 10−2	3.41 × 10−2	1.74 × 10−2	1.38 × 10−1	3.65 × 10−1	4.44 × 10−1	1.48 × 10−2
	gen	1.94 × 10−1	6.42 × 10−2	9.59 × 10−2	9.22 × 10−2	1.95 × 10−1	6.58 × 10−1	7.57 × 10−1	1.34 × 10−1
	model	1.56 × 10−1	1.98 × 10−2	4.68 × 10−2	2.87 × 10−2	1.60 × 10−1	5.44 × 10−1	6.57 × 10−1	2.83 × 10−2
3	plat	4.77 × 10−8	2.11 × 10−3	2.72 × 10−5	1.94 × 10−5	4.37 × 10−3	3.13 × 10−8	1.14 × 10−6	1.20 × 10−2
	gen	9.74 × 10−2	1.62 × 10−2	1.82 × 10−2	1.44 × 10−3	8.00 × 10−3	1.09 × 10−1	3.80 × 10−3	6.24 × 10−3
	model	2.69 × 10−7	1.17 × 10−3	4.09 × 10−5	9.35 × 10−6	1.26 × 10−3	1.86 × 10−7	1.57 × 10−6	2.11 × 10−3

,

Table A7. p-values obtained from two way ANOVA analyses of Red-Edge and water content canopy hyperspectral indices. The values lower than 0.05 have been reported in bold.

cmp	Fact.	RE1	RE2	RE11	RE21	RIRE	NDVI1	REIP1	REIP2	NDWI	NDWI2
1	plat	3.40 × 10−8	5.54 × 10−8	2.34 × 10−7	6.38 × 10−9	1.53 × 10−3	6.95 × 10−7	2.15 × 10−9	5.17 × 10−9	1.09 × 10−2	7.66 × 10−4
	gen	3.60 × 10−1	2.81 × 10−1	1.76 × 10−1	1.44 × 10−1	5.65 × 10−2	3.37 × 10−2	9.94 × 10−2	1.12 × 10−1	6.13 × 10−1	5.03 × 10−1
	model	2.65 × 10−7	4.00 × 10−7	1.38 × 10−6	4.42 × 10−8	1.92 × 10−3	2.26 × 10−6	1.46 × 10−8	3.45 × 10−8	3.87 × 10−2	3.44 × 10−3
2	plat	1.68 × 10−2	3.93 × 10−2	9.40 × 10−2	3.24 × 10−2	5.48 × 10−1	5.09 × 10−2	1.85 × 10−2	1.74 × 10−2	5.21 × 10−1	9.15 × 10−2
	gen	5.21 × 10−2	2.45 × 10−2	2.04 × 10−1	1.08 × 10−1	7.95 × 10−2	6.35 × 10−2	5.04 × 10−2	5.12 × 10−2	3.35 × 10−2	2.49 × 10−2
	model	2.29 × 10−2	2.58 × 10−2	1.28 × 10−1	4.75 × 10−2	1.64 × 10−1	4.98 × 10−2	2.41 × 10−2	2.33 × 10−2	7.76 × 10−2	3.70 × 10−2
3	plat	1.86 × 10−2	1.40 × 10−2	3.89 × 10−8	4.31 × 10−8	5.75 × 10−4	1.19 × 10−2	3.16 × 10−3	3.77 × 10−3	6.42 × 10−6	7.75 × 10−4
	gen	1.33 × 10−2	3.86 × 10−2	1.04 × 10−1	1.24 × 10−1	9.14 × 10−2	9.27 × 10−2	7.44 × 10−1	7.62 × 10−1	3.89 × 10−1	3.93 × 10−2
	model	4.86 × 10−3	8.29 × 10−3	2.26 × 10−7	2.60 × 10−7	1.12 × 10−3	1.30 × 10−2	1.48 × 10−2	1.75 × 10−2	3.89 × 10−5	9.08 × 10−4

,

Table A8. p-values obtained from two way ANOVA analyses of Landsat 8 OLI spectral indices. The values lower than 0.05 have been reported in bold.

cmp	Factor	NDVIL	SAVIL	OSAVIL	EVIS	SIPI	NPCI	NWDI	NWI1	NWI2
1	plat	4.21 × 10−6	1.68 × 10−7	1.01 × 10−6	2.31 × 10−7	3.36 × 10−7	1.86 × 10−5	2.49 × 10−4	5.29 × 10−4	1.01 × 10−6
	gen	1.57 × 10−1	2.58 × 10−1	2.00 × 10−1	2.77 × 10−1	1.35 × 10−1	3.34 × 10−1	2.87 × 10−1	5.68 × 10−1	3.85 × 10−1
	model	1.93 × 10−5	1.12 × 10−6	5.60 × 10−6	1.55 × 10−6	1.80 × 10−6	9.83 × 10−5	9.68 × 10−4	2.62 × 10−3	6.79 × 10−6
2	plat	1.10 × 10−1	1.32 × 10−2	3.71 × 10−2	9.31 × 10−3	1.24 × 10−1	3.59 × 10−1	5.20 × 10−2	8.65 × 10−2	2.25 × 10−1
	gen	2.02 × 10−1	6.19 × 10−2	9.09 × 10−2	5.19 × 10−2	1.81 × 10−1	7.61 × 10−1	1.75 × 10−1	2.72 × 10−2	4.16 × 10−1
	model	1.41 × 10−1	2.06 × 10−2	4.85 × 10−2	1.51 × 10−2	1.44 × 10−1	5.76 × 10−1	7.95 × 10−2	3.85 × 10−2	3.18 × 10−1
3	plat	4.14 × 10−8	6.96 × 10−3	1.57 × 10−4	2.46 × 10−2	2.97 × 10−8	3.77 × 10−6	5.29 × 10−5	4.21 × 10−4	2.42 × 10−3
	gen	8.04 × 10−2	3.20 × 10−1	1.88 × 10−1	4.21 × 10−1	8.48 × 10−2	3.09 × 10−2	1.30 × 10−2	3.61 × 10−2	1.50 × 10−1
	model	2.25 × 10−7	1.81 × 10−2	5.29 × 10−4	6.03 × 10−2	1.67 × 10−7	9.64 × 10−6	5.99 × 10−5	5.33 × 10−4	4.87 × 10−3

,

Table A9. p-values obtained from two way ANOVA analyses of Sentinel 2 MSI structural/LAI and Chlorophyll/pigments spectral indices. The values lower than 0.05 have been reported in bold.

cmp	Factor	NDVI1	NDVI2	NDVI3	SAVI	OSAVI	EVI	SIPI	NPCI
1	plat	5.07 × 10−7	6.29 × 10−7	1.03 × 10−6	2.55 × 10−6	4.01 × 10−6	3.22 × 10−6	6.29 × 10−7	4.39 × 10−9
	gen	1.22 × 10−1	1.38 × 10−1	1.33 × 10−1	2.61 × 10−1	1.98 × 10−1	3.02 × 10−1	1.38 × 10−1	1.08 × 10−1
	model	2.55 × 10−6	3.23 × 10−6	5.05 × 10−6	1.44 × 10−5	1.99 × 10−5	1.87 × 10−5	3.23 × 10−6	2.93 × 10−8
2	plat	1.25 × 10−1	1.16 × 10−1	1.51 × 10−1	1.29 × 10−2	3.42 × 10−2	1.02 × 10−2	1.28 × 10−1	4.08 × 10−1
	gen	1.98 × 10−1	1.92 × 10−1	2.51 × 10−1	6.49 × 10−2	9.30 × 10−2	6.97 × 10−2	1.77 × 10−1	7.99 × 10−1
	model	1.51 × 10−1	1.42 × 10−1	1.92 × 10−1	2.06 × 10−2	4.64 × 10−2	1.77 × 10−2	1.44 × 10−1	6.39 × 10−1
3	plat	3.60 × 10−8	5.11 × 10−7	4.39 × 10−8	6.07 × 10−3	1.35 × 10−4	1.99 × 10−2	5.11 × 10−7	2.13 × 10−6
	gen	6.51 × 10−2	2.93 × 10−1	7.12 × 10−2	3.15 × 10−1	1.80 × 10−1	4.20 × 10−1	2.93 × 10−1	1.24 × 10−2
	model	1.87 × 10−7	3.32 × 10−6	2.30 × 10−7	1.60 × 10−2	4.54 × 10−4	5.06 × 10−2	3.32 × 10−6	4.17 × 10−6

and

Table A10. p-values obtained from two way ANOVA analyses of Sentinel 2 MSI Red-Edge and water content spectral indices. The values lower than 0.05 have been reported in bold.

cmp	Factor	RE11	RE21	NWDI	NWI1	NWI2	NDVI4	REIP1
1	plat	1.58 × 10−5	7.17 × 10−2	7.84 × 10−7	9.53 × 10−4	1.25 × 10−6	2.29 × 10−7	1.00 × 10−6
	gen	8.06 × 10−2	2.27 × 10−1	1.41 × 10−1	4.74 × 10−1	4.54 × 10−1	1.11 × 10−1	4.09 × 10−1
	model	4.88 × 10−5	9.29 × 10−2	3.98 × 10−6	4.06 × 10−3	8.78 × 10−6	1.19 × 10−6	6.88 × 10−6
2	plat	8.54 × 10−2	7.10 × 10−2	4.35 × 10−2	8.79 × 10−2	2.10 × 10−1	1.24 × 10−1	1.88 × 10−2
	gen	1.95 × 10−1	1.57 × 10−1	1.72 × 10−1	2.70 × 10−2	4.13 × 10−1	2.02 × 10−1	5.19 × 10−2
	model	1.17 × 10−1	9.43 × 10−2	6.94 × 10−2	3.87 × 10−2	3.03 × 10−1	1.52 × 10−1	2.47 × 10−2
3	plat	4.14 × 10−8	1.09 × 10−7	6.74 × 10−7	3.69 × 10−4	3.36 × 10−3	2.14 × 10−8	4.33 × 10−3
	gen	1.08 × 10−1	1.09 × 10−1	3.07 × 10−1	3.21 × 10−2	1.45 × 10−1	6.59 × 10−2	7.60 × 10−1
	model	2.42 × 10−7	6.00 × 10−7	4.36 × 10−6	4.50 × 10−4	6.25 × 10−3	1.16 × 10−7	1.97 × 10−2

). The number of discriminated combination pairs of each spectral index in the different measurement campaigns was then used as a proxy of its effectiveness in capturing the effects of the three water salinization levels on the three cardoon genotypes.

Globally, the canopy indices were found to perform better than those assessed by leaf measurements, especially in the most advanced development stages (cmp3), where only the NDVI, NPCI, SIPI RIRE, NDVI1, NDWI1 indices remain sensible to different salinization concentration at leaf level. The behavior of most leaf indices was significantly sensitive to the salinization level (plat label) in the first plant development stages (cmp1, cmp2) while once the plants were sufficiently grown, their capability to capture the effects of different water salinization levels (cmp3) decreases. In general, the leaf spectral indices capability to discriminate the spectral effects of the water salinization level on different genotype is poor, only the NDVI, TVI, SIPI, RIRE, NDVI1, NDWI1 and NDWI2 leaf indices were sensible to genotype effects at leaf level in the first growth stage (cmp1). Their capability to discriminate the effects due to the different cardoon genotypes becomes null (Table 4 leaf reflectance side) in the subsequent development phases (i.e., cmp2, cmp3). Among the foliar spectral indices sensible only to the water salinization effects on cardoon, the most performant, in the 3rd development phase (cmp3) were neatly the RIRE and NDWI (able to discriminate, respectively, two and three of index combinations pairs), with the others at the same level of performance (Table 6). Leaf indices seem mostly ineffective to capture the most subtle responses related to genotype differentiation in more advanced development stages coupled with a poorer effectiveness in detecting spectral effects of water salinization. In fact, at leaf level, only the RIRE was found capable of discriminating the spectral effects of all three levels of water salinization corresponding to three growing plots in the 3th measurement campaign (Table 4).

Contrary to their leaf counterparts, the canopy spectral indices discrimination power of spectral effects linked, respectively, to water salinization level and genotype, increases with the development stage, according to the amplification of the spectral effects of the 3-d vegetation architecture growth compared to those of single leaves. In fact, all the models related to the selected canopy spectral indices became significant in the last development stage (cmp3, canopy reflectance side), with the effective ability of all indices to discriminate the water salinization effects, coupled with higher ability for 9 of them to capture also spectral effect differentiation related to the three cardoon genotypes (Table 4, cmp3). The performance of these latter indices, evaluated for the cmp3 data by Tukey’s test, was reported in the canopy reflectance side of Table 6. These achievements demonstrate that the OSAVI and TCARI indices were able to discriminate the spectral effects related to all water salinization levels, while SAVI, TCARI, and TVI were the most sensible to differentiations in the genotype spectra, with two pairs discriminated. Finally, the most effective index in discriminating the different levels of salinization of the water and the different cardoon genotypes was the TCARI, with its total score of five (three water salinization level + two genotype combination pairs discriminated). In the “total” column of Table 4, the total sum of the combination pairs discriminated was reported as proxy of global performance of the index (for both genotype and plateau factors at leaf and canopy level). Between the structural indices, the SAVI scored the maximum number of discriminated pairs for both plateau and genotype factors, at leaf and canopy level, while TCARI resulted in the best performing of chlorophyll/pigment indices with a superior score number. Globally, also the red-edge indices have shown a good effectiveness with RE2 and REIP2 that maintained a significant performance in discriminating spectral effects of water salinization for different genotypes in the second measurement campaign. Despite its scarce capability to detect the spectral effects due to different genotypes, the NDWI of water spectral indices resulted in the best performing respect to all others. The global score, as number of discriminated pairs for both plateau and genotype combination pairs, at leaf and canopy levels, assessed for each spectral index ranges between 7 and 13, with many indices of the four groups gathering an intermediate value (9–10).

Most of L8 OLI structural indices have shown their effective capability to detect the spectral differences due to effects of different water salinization levels during all the plant growth phases (Table 5). The other chlorophyll/pigment and water indices performed well only in the early and late stages. Their discrimination capability of spectral differences due to different genotypes is weak and limited to the advanced development stage. The SAVI_L, EVI_L, and NPCI_L were the best overall in terms of total number of discriminated pairs.

The Sentinel 2 MSI structural indices showed mostly a similar trend in detecting water salinization spectral effects (Table 6) to which we added the useful capability of those based on Red-Edge Inflection Point (REIP) working also in intermediate and advanced stages of plantar development. The capability to detect the most subtle spectral differences in the three cardoon genotypes due to water salinization levels (gen) remains weak and limited to intermediate and advanced plant growth stages. The most performant indices, which reached a total score of six combination pairs discriminated, were SAVI_S, NPCI_S, NWI1_S, and the red-edge-based that demonstrated also the ability to detect genotype spectral differences both at the intermediate and advanced stages of plant development.

The different leaf and canopy hyperspectral indices through the data acquired during the three measurement campaigns evidenced a wide range of spectral indices sensibility, at plateau and genotype levels. In the early stage (cmp1), both the leaf and canopy indices demonstrated a useful capability to detect the plant spectral differences induced by the three water salinization levels, while only water coupled with some structural and pigment indices was found able to capture the effects at the genotype level. In the following stage (cmp2), contrary to canopy indices, the leaf ones maintained their aptitude to discriminate the effects of water salinization level without genotype distinction.

In the later stage (cmp3), the sensibility of the canopy indices, especially those of the structural and pigment ones including the genotype effects increased, while diminishing the reactivity of the leaf indices decreased. Similar trends were found for the S2 MSI broadband indices with additional capabilities demonstrated by those based on the red-edge, acquisition channels (Figure 4). Both the OLI and S2 broadband indices demonstrated a suitable aptitude to detect the spectral effects of water salinization levels on cardoon plants, especially in early and sufficiently developed stages (1.2–1.4 stages of BBCH scale), with a weak additional sensitivity to different genotypic effects of water indices in these latter.

3.3. Biometry Modelling

Figure 6 shows the boxplot graphs related to Biomass (mg/cm²) and LAI (m²/m²) biometric measurements derived from cardoon samples acquired on field during the second campaign (cmp2). The graphs a and b show data related to the three plateaus with different water salinization levels (starting from A with lowest salinization to C with the maximum one). The graphs c and d show the boxplots assessed for the three different C. cardunculus genotypes.

The average trend of LAI and biomass density (Bio) do not show significant (at 95% of confidence level) variations between plateaus A and B, with raising water salinization, while the differences get significant for samples of C plateau (p-value < 0.05). This find is confirmed also by the results of Tukey’s test in

Table 7. Two-tails ANOVA and Tukey’s test analysis of biometry data. In bold have been reported the entries significant at 95% of confidence level (p-value < 0.05). The plat and gen labels indicate, respectively, the plateau (different water salinization) and cardoon genotype factors, while models refer to the ANOVA global model. The Tukey entries refer to three combination couples of plateau and genotypes).

Parameter	p-Value			Tukey’s Test
				Plat				Gen
	Plat	Gen	Model	A-B	A-C	C-B	n. < 0.05	g1-2	g1-3	g2-3	n. < 0.05
LAI	5.0457 × 10−8	0.55587	4.2705 × 10−7	0.11605	1.8986 × 10−8	4.9028 × 10−6	2	0.94405	0.73846	0.54163
Bio	0.00653	0.35543	0.01826	0.97825	0.01135	0.01785	2	0.99712	0.44424	0.40467

, where only the two pairs A-C and B-C of LAI/Bio mean values were discriminated at plateau level, while no significant variation was found for samples of different genotypes (p-value > 0.05).

In agreement with these results, most of canopy indices of the second campaign show significant differences only at plateau level, while a weak discrimination of genotype effects was shown only by red-edge and water spectral indices (Table 6, canopy reflectance, cmp2). Similar trends were found for the OLI and MSI broadband indices.

Figure 7 shows the plot of the biomass density (Bio) and the corresponding LAI measurements including their best-fit allometric curve as:

$$\mathrm{LAI}=6.4141{\mathrm{Bio}}^{-0.954}$$

(3)

with the significant correlation found between them evidenced by an R²_adj (adjusted correlation coefficient) = 0.53197. The LAI shows a rising trend linked to that of the lowering Bio, depicted by the fitting allometric curve, with the best-fit intermediate LAI values around 1.25. The Bio variable measurements registered the maximum in case of the smallest leaf development caused by water salinization increase, in particular for the 3rd genotype whose values are the highest ones. The distribution of measurements from the samples of the C plateau (triangular dots) is mainly located in the area near the origin of the axes, in the lower left corner of the graph (average LAI = 0.112). The other data, referring to decreasing water salinization levels, spread along the best-fit curve, with those related to the B plateau (mean LAI = 1.175) and those from the A plateau above them (mean LAI = 1.518). From this graph and from the previous analyses of the biometric data, it appears that, above a water salinization threshold (i.e., water salt concentration of the plateau C), the global development of cardoon plants is significantly inhibited, while under this threshold, the water salinization neither affects the development of the leaf surface (LAI) nor the biomass density (Bio) of C. cardunculus plants. In order to link the biometric data to the corresponding spectral indices for modelling purposes, both the average values at genotype and plateau levels have been assessed from the related measurement datasets. In particular, the leaf and canopy hyperspectral indices and the broadband ones of the second campaign have been selected to be extensively tested for linear modelling of the yield Y, derived from measured LAI and Bio biophysical parameters (Equation (1)). The results of this modelling approach were reported below, under the form of bi-dimensional scatter plots where the coordinates of each index are the R correlation (Spearman) coefficients derived from the linear model of Y (y-axis) and LAI (x-axis) parameters, more directly linked to biomass resources on field by the bioenergetic crops. In particular, Figure 8a shows the scatter plots for canopy hyperspectral indices, while Figure 8b includes that assessed for the Sentinel MSI broadband indices. In the graphs, some spectral indices that overlap each other are not reported to avoid confusion and improve the clearness of distributions. In general, in all the graphs, the position related to all indices is distributed along the bisector of the second and third quadrants, which means that the correlations of the spectral indices with LAI and Y may have the opposite sign depending also on the implicit contribution of Bio in Y without significant differences among the three genotypes.

As reported in graph a of Figure 8, the correlation of the red-edge (REIP, RE21, RE11) canopy’s spectral indices, with the LAI/Y is positive and significant, with that slightly inferior of NDVI, NDVI1, and SIPI. The soil corrected SAVI and OSAVI, RIRE, and the NDWI indices showed instead poor correlation values with those of SAVI inferior of others and negative for both LAI and Y. The other indices (MCARI and NPCI) are characterized by a varying negative correlation with TCARI reaching the minimum values.

Most of the S2 broadband spectral indices (Figure 8, graph b), presented a behavior similar to that of the corresponding hyperspectral canopy index, with those red-edge-based (RE11, RE12, REIP) more concentrated and having the maximum positive correlation with LAI and Y. Contrary to NDVI, the EVI index, specifically designed for HR satellite optical sensors (to minimize soil and atmospheric noise contribution), shows a negative correlation, probably due to the overcorrection of absent atmospheric effects, based on blue channel. The correlation of the OLI spectral indices showed a trend similar to that of the S2 MSI, with the NDVI achieving the maximum positive correlation but inferior to that of S2 red-edge ones. In general, from these preliminary R correlation test, the red-edge indices evaluated for the second field campaign (cmp2) happened most suitably in modelling the Y biometric data (Figure 8a,b), maintaining at same time a discrimination capability for the corresponding effects of different water salinization levels (Table 6, canopy reflectance indices, cmp2).

The most relevant positive correlation with Y and LAI is provided by the RE21 index, normalized difference red-edge index that reached the maximum value, followed by RE11 (Figure 8). The biomass production assessment at field level requires spatial modeling of the yield Y, in order to provide input for sustainable management and bioenergy exploitation. In this perspective, a deeper analysis of this index was carried out considering that currently the hyperspectral and multispectral data collected remotely of the areas of interest for spatial distributions assessment can be provided by the airborne or satellite operating platforms (i.e., S2 MSI satellite). Figure 9 shows the average trends of the Y derived from the biometric data and the corresponding RE21 values for each plateau (different dot shape) and genotype (different dot color). The regression equation found is:

$$\mathrm{Y}=55.15033\mathrm{RE}21-7.98479$$

(4)

The R² correlation of the related Y-RE21 linear model scored 0.836 (R²a_dj = 0.656) that means an effective capability to capture the above-ground biomass spatial variations supported by a model significance better than 95% (F-value = 16.26, Prob > F = 0.00498). The lowest Y (around 1.5 t/ha) and RE21 values were found for the cardoon plants of the three genotypes of the C plateau (triangular dots) fed with the most salinized water. The maximum (around 9 t/ha), instead, has been reached both by the first cardoon genotype grown in water at intermediate salinization (circle black dot) or not salinized water (square black dot) and by the second genotype (circle red dot) irrigated with water at intermediate salt concentration. The Y of plants of genotype 1 (black dots) and 2 (red dots), grown using water at intermediate salt concentration, approaches 9 t/ha (maximum); in particular, the Y of genotype 2 exceed that obtained with water non-salinized (about 8 t/ha).

The RE21 index and the corresponding Y values of all cardoon genotypes feeding with the most salinized water are rightly the lowest (triangular dots), but they differ in the trends of the remaining plateau. In fact, although the Y rising trend of the genotype 1 (black dots) is in agreement with water salinization reduction, the Y and RE21 index values corresponding to water middle salinization level (circle dot) were found too high respect to those of the genotype 3, while the genotype 2 showed a similar Y value compared to a significantly lower index. All the plateau averages of the genotype 3 (green dots) are rightly located along the fitting line, with Y related to middle water salinization significantly smaller than the other two and that obtained through the non-salinized water located between the others corresponding two. All the plateau Y averages of the second genotype (red dots) are located above and at a significant distance from the red line of best fit with the yield Y corresponding to the intermediate water salinization (red circle dot) higher than that of the plants fed with no salinized water (red square dot), which present an inversion not evidenced by the trends of the measurement of the other two genotypes.

A similar trend was found for the RE21_s index, derived from the S2 MSI broadband spectral responses, with a higher R² correlation coefficient and the remaining modelling parameters reported in

Table 8. Y linear model statistical parameters using multispectral S2 MSI RE21 index.

-	DF	Sum of Squares	Mean Square	F Value	Prob > F
Model	1	132.77595	132.77595	14.98502	0.00612
Error	7	62.02406	8.86058
Total	8	194.80001

, including those stating the relevant model confidence level.

Although the ANOVA analysis of biometry data did not evidence significant differences of genotypes’ responses to three water salinization levels, the Y biomass growing trends of genotypes 1 and 3 were found in raw agreement each other and with water salinization decreasing, except for intermediate values, corresponding to different rises in Y and corresponding RE21 index. In particular, the Y and RE21 intermediate values of this genotype (black circle dot) approaches those referring to those obtained for no salinized water (black and green square dot). The plants of 2nd genotype (red dots) instead, after a more relevant growth of the Y related to intermediate water salinization (similar to that of genotype 1 but without a corresponding rise in RE21 index) show a decreased Y respect to the same genotype thistle grown with water not salinized. In addition, these plants presented a rise in corresponding RE21 (square red dot), equivalent to around one half respect to others without water salt concentration. These results show that, at irrigation water salt concentrations of 200 mM/L (plateau C), the growth of all thistle genotypes used here is significantly inhibited both in term of Y and spectral indices responses (Y = 1–2 t/ha), while at half concentrations (plateau B), the different genotypes exhibit different behaviors. In particular, the plants of genotype 1 and 2 (black and red circle dots) evidence productivity growth (Y = 9 t/ha) approaching those obtained for genotype 1 and 3 (black and green square dots), with water not salinized (plateau A). The Y values corresponding to the plateau A for genotypes 1 and 3 look quite similar, while that of the 2nd genotype diminishes respect the corresponding Y amount of plateau B (100 mM/L).

Ultimately, from these preliminary findings, the intermediate salt concentration around 100 mM/L seems to have a negative impact on the productivity Y of the cardoon genotype 3, not affecting that of the genotype 1 and favoring the development of the second genotype used here. Although the intermediate water salinization feeding involves a rise in RE21 spectral index responses of plants of all three genotypes, their growth amount differs. The plants of the cardoon genotype 2 exhibit the smallest index value, those of the cardoon genotype 3 presented intermediate RE21 values, while the highest were assessed for the genotype 1 plants. From the comparison with the data referred to samples irrigated without salinized water, we can evidence that the 100 mM/L concentration does not impact on the Y of cardoon genotype 1, while it improves productivity of the genotype 2 (of about 1 t/ha) and inhibits that of the genotype 3 (of about 3 t/ha), as can be derived from the differences in y-coordinates in the graphs. In any case, it should be highlighted that, although Y (t/ha) is specifically linked to the production of extensive field cultivation, its estimated maxima values in our controlled and limited experimental context (about 10 t/ha of above-ground biomass) are roughly compatible with those reported in the literature for bioenergy productions growing in the field at comparable agronomic condition and with low fertilizer inputs, taking into account also of biomass partitioning between leaves, stalks, and heads and wet/dry mass mutual dependence [6,9].

4. Conclusions

Cardoon exploitation as a bioenergy resource is widely diffusing, particularly in the Mediterranean basin, due to its relevant annual biomass yield with a low request of inputs and a significant resilience to stresses due to drought and water shortage, soil/water salinization, and soil pollution. This crop may have a wide range of bio-industrial applications and fit into different intensive cultivation strategies that can be improved by precision farming approaches, supported by the currently available remote and proximal sensing techniques, including those exploiting EO data continuously provided by orbiting satellites. In this perspective, various spectral vegetation indices based on these remote sensing techniques tested here have demonstrated their ability to discriminate between the spectral differences induced on cardoon plants by different water salinization levels (at leaf and canopy level) at the different growth stages. The more difficult detection of subtler spectral effects of water salinization on the different cardoon genotypes was anywhere provided by different indices in the various development phases. Even if, as expected, the tested broadband spectral vegetation indices showed a slightly lower capability, they maintained an appreciable discrimination of the spectral effects due to three salinization levels in irrigation water of cardoon growing plants with a weak capability to differentiate the genotype effects. The red-edge-based indices demonstrated also a good modelling capability of biometric parameters derived from the field data as the yield, Y (t/ha). The useful cardoon Y (t/ha) modelling capability, demonstrated by narrow/spectral vegetation indices broadband tested here (in particular by red-edge-based ones), in perspective, can support effective, repetitive, and extensive monitoring of these bioenergetic crops by providing early detection and mapping of plan stresses and health that can have an impact on production. In the framework of precision farming approach, these information are the basis for sustainable and effective management with space-variant interventions (even supported by semi-autonomous agricultural VRT—variable rate technology machinery) for optimizing the production minimizing at same time external inputs and environmental impacts. The current wide availability of HR satellite multispectral data, remotely detected in the various spectral ranges, including the red edge one (i.e., S2 MSI), with the realistic perspective of exploiting the available satellite hyperspectral sensors (PRISMA, ENMAP), provide very useful monitoring opportunities of cardoon and other bioenergy crops. Although the noise effects due to atmospheric disturbance and the sun-view anisotropy (minimized through acquisition geometry and normalized spectral indices exploitation) have not been taken into account here, the satellite broadband indices in different wavelengths evaluated have preliminarily proven to be potentially useful for extensive and effective cardoon monitoring. In addition, the remote/proximal sensing techniques tested here could support the further extensive characterization of cardoon production, including the quantification of seeds and biochemical substances and compound concentrations within a biorefinery context, with the selection of suitable genotypes for biomass production also characterized by a different degree of the chlorine absorption, which is a further undesirable factor to be taken into consideration for the processes for bioenergy production.

5. Reference Projects

The present research has been conducted within the framework of:

• Multi-year national program of National Electric System Research (RdS) financed by the MISE (Ministry of Economic Development) for research and development activities in the electricity sector and to improve its cost-effectiveness, safety, and environmental compatibility, with funding to ENEA for the development of models and technologies to assess and reduce anthropogenic impacts and natural hazards arising from energy production, in particular for the implementation of innovative monitoring methods based on Earth observation techniques;
• REBIOCHEM project, funded by the Italian Ministry of Education, University and Research for the cluster on Green Chemistry, focusing on assessing the potential of different cardoon species and genotypes in terms of biomass production and extraction of useful compounds from raw materials, within the biorefinery context.