Hyperspectral and Thermal Sensing of Stomatal Conductance, Transpiration, and Photosynthesis for Soybean and Maize under Drought

: During water stress, crops undertake adjustments in functional, structural, and biochemical traits. Hyperspectral data and machine learning techniques (PLS-R) can be used to assess water stress responses in plant physiology. In this study, we investigated the potential of hyperspectral optical (VNIR) measurements supplemented with thermal remote sensing and canopy height (h c ) to detect changes in leaf physiology of soybean (C 3 ) and maize (C 4 ) plants under three levels of soil moisture in controlled environmental conditions. We measured canopy evapotranspiration (ET), leaf transpiration (T r ), leaf stomatal conductance (g s ), leaf photosynthesis (A), leaf chlorophyll content and morphological properties (h c and LAI), as well as vegetation cover reﬂectance and radiometric temperature (T L,Rad ). Our results showed that water stress caused signiﬁcant ET decreases in both crops. This reduction was linked to tighter stomatal control for soybean plants, whereas LAI changes were the primary control on maize ET. Spectral vegetation indices (VIs) and T L,Rad were able to track these di ﬀ erent responses to drought, but only after controlling for confounding changes in phenology. PLS-R modeling of g s , T r , and A using hyperspectral data was more accurate when pooling data from both crops together rather than individually. Nonetheless, separated PLS-R crop models are useful to identify the most relevant variables in each crop such as T L,Rad for soybean and h c for maize under our experimental conditions. Interestingly, the most important spectral bands sensitive to drought, derived from PLS-R analysis, were not exactly centered at the same wavelengths of the studied VIs sensitive to drought, highlighting the beneﬁt of having contiguous narrow spectral bands to predict leaf physiology and suggesting di ﬀ erent wavelength combinations based on crop type. Our results are only a ﬁrst but a promising step towards larger scale remote sensing applications (e.g., airborne and satellite). PLS-R estimates of leaf physiology could help to parameterize canopy level GPP or ET models and to identify di ﬀ erent photosynthetic paths or key physiological variables regulating water and carbon assimilation: leaf stomatal leaf and leaf M.G., T.N.M., and A.B.; Camera calibration, M.G., C.J.K., L.G., and V.S.-P.; Statistics: M.G., A.B., and V.S.-P., Resources, M.S.J., T.N.M., and M.G.; Results analysis, V.S.-P., M.G., T.N.M., A.B., M.S.J., S.L., and X.M.; Supervision, M.G., T.N.M., S.L., and X.M.; Writing—original draft preparation, V.S.-P.; Writing—review and editing, all authors. Visualization, V.S.-P. All authors have read and agreed published of the manuscript.


Introduction
Irrigation for agriculture represents the most significant demand for freshwater uptake of around 70% of all water withdrawal [1,2]. This demand is growing due to climate change. Increases in the intensity and frequency of droughts and heatwaves will affect soil water demand and supply of crops, hindering crop production [3][4][5]. Crops such as maize and soybean, sensitive to droughts, are already experiencing reductions in yields [6][7][8]. These staple crops are extensively grown for food production worldwide, covering areas of 139 million ha for maize and 75.5 million ha for soybean crops [9]. Thus, further information of their responses to droughts at the farm level can improve the sustainability of water use.
In response to soil drought, or water stress, plants react with various physiological, biochemical, and morphological adjustments, observable early at the drought onset (fast changes) or after some time (slow changes) [10,11]. Generally, as a result of soil drought, leaf water potential (Ψ L ) decreases and stomata close, which limits transpiration (T r ) and photosynthesis (A) to different degrees depending on plant type, life history and environmental factors [12][13][14]. Thus, stomatal conductance (g s ) is often used as an indirect proxy to detect plant water status and stress [12,15,16]. Cell enlargement is also a very sensitive process to mild water stress [11,13], which it will be eventually reflected in lower canopy heights (h c ) or leaf area index (LAI) [7,17]. In addition, there are other delayed responses in pigments such as reductions in chlorophyll content [7,10,13,18] or carotenoids [19]. Chlorophyll reduction caused by water stress can be attributed to chlorophyll degradation, pigment photo-oxidation, and insufficient synthesis of chlorophyll [17]. However, in the assessment of crop physiological responses to drought, it is not always easy to separate drought effects from phenology and leaf age effects [20]. Accounting explicitly for different hydraulic and metabolic traits is necessary to improve predictions to soil water stress. For example, maize plants have C 4 photosynthetic pathway, while soybean are classified as C 3 , being less efficient under warm temperatures as their photorespiration is higher [21,22]. Soil water stress could lead to more substantial declines of g s in C 3 species than C 4 , due to metabolic limitations [23,24]. In addition, species with more anisohydric behavior, as reported for maize [25] can delay stomatal closure, maintaining photosynthesis level. In contrast, more isohydric plants will rapidly close stomata in response to water deficit, limiting water stress damage but reducing CO 2 assimilation [16]. Incorporation of some of these different traits or responses in models (e.g., Priestley-Taylor Jet Propulsion Laboratory (PT-JPL) [26][27][28], Light Use Efficiency Gross Primary Productivity (LUE GPP) [26,29,30] and Soil Canopy Observation, Photochemistry and Energy fluxes (SCOPE) [31]) is key for better predictions and crop management.
Remote sensing techniques have been widely used as a non-destructive approach to estimate plant biochemical and morphological constituents that modify energy absorption and scattering in different spectral wavelengths [20,32,33]. Relevant traits, such as hydraulic safety or photosynthetic pathways, can be detected through optical or thermal features [10,34]. However, the use of remote sensing to estimate plant function or physiology is very complex, as the mechanisms linking reflectance and emission with plant functional gradients are not always explicit or known [35]. In remote sensing applications, the most frequent method to monitor plant responses to soil water stress is using vegetation indices (VIs) [10,15,[35][36][37]. VIs can be associated to plant structural properties like LAI (e.g., normalized difference vegetation index (NDVI) [38]) or pigment changes such as chlorophyll (e.g., transformed chlorophyll absorption in reflectance index (TCARI) [39]) and carotenoids (e.g., pigment specific normalized difference index (PSNDc) [40]). The photochemical reflectance index (PRI) is a more direct estimate of plant physiology, tracking photosynthetic efficiency [41], but it is also affected by changes in LAI or viewing geometry that might obscure the effects of drought [42,43]. Because leaf temperature usually increases rapidly under drought due to less transpiration cooling, thermal imaging has been applied to monitor drought responses [37,[44][45][46][47]. Nonetheless, canopy temperature is also influenced by radiation budgets, atmospheric conditions, and leaf and canopy traits [13,48], which makes it difficult to isolate the effect of soil water stress. Growth chambers provide the possibility to set and control some of these factors such as air temperature, humidity, radiation, wind speed, and soil moisture [49]. For example, thermal variation was observed between water stress and well-watered maize plants at the seedling stage in a growth chamber [50], showing canopy temperature as a suitable and leading indicators to detect water stress.
Combination of different spectral regions provides spatial and temporal data on crop water use [10], which is essential to simulate crop productivity [46,51]. Synergies between different remote sensing domains (e.g., visible (VIS), near-infrared (NIR), fluorescence, thermal, and microwave) have been proven to better model and estimate water use efficiency (WUE), evapotranspiration (ET), and gross primary productivity (GPP) [52]; detect water stress patterns [37,[45][46][47]; and monitor crop yield [53]. Hyperspectral imaging sensors with contiguous narrow bands along the visible and near infrared (VNIR) domain allow to apply spectroscopy and chemometrics techniques at high spatial resolution, improving the accuracy and types of physiological variables retrieved [15,18,20,54]. In addition, it explores new relations between spectral variables and plant function [34]. Machine learning methods such as the partial least squares regression (PLS-R) can handle multiple and correlated spectral bands, being more suitable for large datasets outperforming physiological estimates from spectral indices [15,20,54,55]. Furthermore, from PLS-R models, we can identify the most relevant wavelengths contributing to the model, which might not be necessarily those from commonly used VIs [15]. PLS-R has been successfully used to predict physiological properties such as A [56,57], g s [15,57], T r [57], maximum Rubisco activity (V cmax ) [58][59][60], and Ψ L [15]. PLS-R has been also applied to model biochemical parameters such as chlorophyll content [56,61], carotenoids [56], and nitrogen concentration [56,61,62], and morphological parameters (e.g., leaf mass per unit area [56,61,62] and crop yield [53]). These studies have been conducted across species (e.g., wheat, vineyards, soybean, maize, pasture, rice, and boreal trees) over Mediterranean to tropical climates and in greenhouse or external field conditions. Although using remote sensing to assess crop responses to drought is a very active topic of research, most studies to date have focused on estimating drought responses looking at biochemical and structural parameters. In this study, we aim to assess physiological responses to soil drought for two crops with different photosynthetic pathways and water-use strategies, using hyperspectral optical and thermal remote sensing under controlled conditions while separating phenology effects during a vegetative period. Specific objectives for this study were to: • Assess responses to soil water stress of soybean and maize physiology (canopy ET, leaf g s , leaf T r and leaf A), morphology (h c , LAI), biochemistry (chl), and remote sensing (T L,Rad and reflectance) variables. • Investigate PLS-R as a method to model leaf g s , T r , and A using hyperspectral data and identify the most relevant wavelengths contributing to the model as well as the model improvement when incorporating T L,Rad and h c into the models.
The results from this study can be useful for applications in growth chambers and outdoors at the farm level using unmanned aerial vehicles (UAVs) carrying hyperspectral, thermal cameras and obtaining h c from structure from motion (SfM) photogrammetry or light detection and ranging (LiDAR) sensors.

Experimental Design
The experiment was carried out between March and June of 2018 in a controlled growth chamber at the Risø Environmental Risk Assessment Facility (RERAF) phytotron (Technical University of Denmark (DTU), Roskilde campus, Denmark), with a ground area of 24 m 2 and height of 3.1 m (e.g., [63,64]). Air temperature in the chamber was set to 25 • C for daytime and 19 • C for night-time. The relative humidity (RH) was kept at 50% during day hours and 70% during night time. To promote rapid growth, illumination conditions simulated 16 hours of daytime. Light racks were composed by 28 high-pressure mercury lamps (1000 W) and 14 halogen lamps (400 W). CO 2 concentration was set constant at 400 ppm, similar to the global average. Air temperature, relative humidity, light intensity, and CO 2 concentration were automatically recorded using a climate control computer.
In the experiment, two types of crop seeds were sown: Zea mays cv. NK Arma (maize) and Glycine max (L.) cv. Merr. Buenos (soybean), provided by SAGEA Centro di Saggio s.r.l. Integrated Solutions for Sustainable Agriculture (https://www.sagea.com). These crops were selected as they are staple crops produced under similar climatic conditions and present different photosynthetic pathways (C 3 for soybean and C 4 for maize) and hydraulic strategies. Soybean seeds were sown on 23 March and maize seeds on 9 April in 9 L pots with 4 kg of soil. Two plants were grown in each pot. Table 1 shows days after sowing (DAS) and days of measurement (DoM) of each crop. Along with this document, DoM was used for visually more precise figures and as reference for phenology. The soil used was a pre-fertilized soil extracted from peat blocks and nutrient enriched (Pindstrup substrate no. 6, Pindstrup Mosebrug A/S, Ryomgaard, Denmark, https://www.pindstrup.com). This soil type has been previously used in other experiments in the same growth chamber [63,64] and presents a total water holding capacity (WHC) of 0.625 kg water/kg soil. Three soil water levels were established to maintain: (i) 100% WHC (control/wet), (ii) 70% WHC (mid), and (iii) 40% WHC (dry), later referred to as soil drought or water stress. Irrigation was applied every three days until the 24 April and after that, every two days. For each WHC group and crop, six replicates were arranged for a total of 36 pots (2 crops × 3 WHC × 6 replicates).

Gas Exchange and Plant Measurements
The gravimetric soil water content was measured by weighing each pot with an electronic balance before and after irrigation. Canopy evapotranspiration (ET, mm day −1 ) was calculated from the water balance to establish the amount of water to be dispensed during each irrigation.
Physiological, morphological, and biochemical measurements were conducted seven times for maize plants and eight times for soybean plants (Table 1) during a time span of about 20 days between 11:00 and 14:00 after irrigation. We measured gas exchange of CO 2 and H 2 O at the leaf level with the open photosynthesis system LI-6400 (Li-COR Biosciences Inc., Lincoln, NE, USA). The system maintained the leaf chamber at the ambient temperature, at 400 µmol mol −1 CO 2 level and constant photosynthetic active radiation (PAR) of 1500 µmol m −2 s −1 (saturation point). Gas exchange measurements were conducted on two randomly selected upper leaves, fully matured and light exposed, to obtain parameters such as photosynthetic CO 2 assimilation rate (A, µmol CO 2 m −2 s −1 ), transpiration rate (T r , mmol H 2 O m −2 s −1 ), stomatal conductance (g s , mol H 2 O m −2 s −1 ), leaf vapor pressure deficit (VPD leaf , KPa), intercellular CO 2 concentration (C i , µmol CO 2 mol −1 moist air), and water use efficiency (WUE = A/Tr, µmol CO 2 mmol −1 H 2 O). We measured absolute chlorophyll concentration (chl, µg cm −2 ) with the SPAD 502 Plus leaf chlorophyll meter (Konica Minolta Sensing Inc., Sakai, Osaka, Japan) in three randomly and fully expanded top leaves. In each leaf, we took the average of three sampling points distributed over the leaf. Selected leaves for chl measurements were not necessary the same as leaves selected for gas exchange measurements. We manually measured canopy height (h c , m) with a meter from the soil surface to the highest point of each plant. h c was used as a proxy for plant growth, as it can be associated to cell enlargement [7,17]. Gas exchange, chl and h c measurements were averaged per pot.
Grouping pots by crop and soil water status, we measured leaf area index (LAI) three times with the plant canopy analyzer LAI-2200C (Li-COR Biosciences Inc., Lincoln, NE, USA) on the last day of the experiment (15 June 2018).

Remote Sensing Measurements
Hyperspectral images were obtained with the Cubert UHD 185 camera (UHD; Cubert GmbH, Ulm, Germany, https://cubert-gmbh.com) with spectral bands between 450 and 950 nm (visible and near infrared ranges, VNIR) along 125 channels with mean spectral resolution of about 10 nm in the VIS, 20 nm in the red-edge, and 40 nm in the NIR region. The camera spatial resolution is 50 × 50 pixels. The focal length of the lenses is 23 mm and the Field of View (FoV) in both directions is 15 • . The Cubert also provides a panchromatic image of 1000 × 1000 pixels.
Thermal images were gathered from a FLIR Tau2 324 (FLIR Systems Inc., Wilsonville, OR, USA), detecting longwave radiation between 7.5 µm to 13.5 µm. The spatial resolution is 324 × 256 pixels with a focal length of 9 mm, which provides a FoV of 48.5 • in the horizontal side and 39.1 • in the vertical view. Figure 1 shows the experiment setup. Both cameras were located approximately 2 meters above the pots and connected to a computer triggered manually after gas exchange measurements, resulting in a pixel size of 10.5 mm/px for the hyperspectral camera and 6.1 mm/px for the thermal sensor. Cameras were positioned in a payload box usually carried by a UAV as shown in Wang et al. (2019) [52]. Due to instrumentation requirements and space management, the separation of cameras was about 18.6 and 3.7 cm in the horizontal axis.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 32 used as a proxy for plant growth, as it can be associated to cell enlargement [7,17]. Gas exchange, chl and hc measurements were averaged per pot. Grouping pots by crop and soil water status, we measured leaf area index (LAI) three times with the plant canopy analyzer LAI-2200C (Li-COR Biosciences Inc., Lincoln, NE, USA) on the last day of the experiment (15 June 2018).

Remote Sensing Measurements
Hyperspectral images were obtained with the Cubert UHD 185 camera (UHD; Cubert GmbH, https://cubert-gmbh.com) with spectral bands between 450 and 950 nm (visible and near infrared ranges, VNIR) along 125 channels with mean spectral resolution of about 10 nm in the VIS, 20 nm in the red-edge, and 40 nm in the NIR region. The camera spatial resolution is 50 × 50 pixels. The focal length of the lenses is 23 mm and the Field of View (FoV) in both directions is 15°. The Cubert also provides a panchromatic image of 1000 × 1000 pixels.
Thermal images were gathered from a FLIR Tau2 324 (FLIR Systems Inc., Wilsonville, OR, USA), detecting longwave radiation between 7.5 µm to 13.5 µm. The spatial resolution is 324 × 256 pixels with a focal length of 9 mm, which provides a FoV of 48.5° in the horizontal side and 39.1° in the vertical view. Figure 1 shows the experiment setup. Both cameras were located approximately 2 meters above the pots and connected to a computer triggered manually after gas exchange measurements, resulting in a pixel size of 10.5 mm/px for the hyperspectral camera and 6.1 mm/px for the thermal sensor. Cameras were positioned in a payload box usually carried by a UAV as shown in Wang et al. (2019) [52]. Due to instrumentation requirements and space management, the separation of cameras was about 18.6 and 3.7 cm in the horizontal axis. Figure 1. Camera set-up and crop overview inside the growth chamber. Two crops were cultivated: soybean and maize, in three water holding capacity (WHC) groups: dry (40% WHC), mid (70% WHC), and wet (100% WHC). For each set we grew six replicates. Images of plants were taken by moving each crop under the box containing thermal and hyperspectral cameras. The wooden box was located 2 m above the pots. Dark blue dotted/dashed line shows the field of view (FoV) of the thermal Flir camera while the red dashed line shows the FoV of the hyperspectral Cubert camera. Both cameras were controlled with an external computer located inside the chamber. Light racks were composed by 28 high-pressure mercury lamps (1000 W) and 14 halogen lamps (400 W). Air temperature and RH were automatically recorded using a climate control computer located outside the growth chamber. Figure 1. Camera set-up and crop overview inside the growth chamber. Two crops were cultivated: soybean and maize, in three water holding capacity (WHC) groups: dry (40% WHC), mid (70% WHC), and wet (100% WHC). For each set we grew six replicates. Images of plants were taken by moving each crop under the box containing thermal and hyperspectral cameras. The wooden box was located 2 m above the pots. Dark blue dotted/dashed line shows the field of view (FoV) of the thermal Flir camera while the red dashed line shows the FoV of the hyperspectral Cubert camera. Both cameras were controlled with an external computer located inside the chamber. Light racks were composed by 28 high-pressure mercury lamps (1000 W) and 14 halogen lamps (400 W). Air temperature and RH were automatically recorded using a climate control computer located outside the growth chamber.

Image Pre-Processing
Hyperspectral Images A radiometric and spectral calibration of the Cubert hyperspectral camera was performed at DTU photonics lab. The spectral response of each pixel was obtained using a monochromator from 430-1030 nm at 1 nm intervals. Subsequently, a radiometric calibration was performed using an integrating sphere to obtain a calibration factor for each pixel in each specific band under nine different illumination levels. The camera integration time was set from 0.5 to 19.9 ms (0.1 ms intervals) at nine different illumination levels. In the growth chamber, we used a spectralon target to characterize the average incoming spectral radiation from the lamps and its fluctuation. On five dates, we acquired a hyperspectral image of the spectralon before measuring the crops. For the remaining dates, we used the average irradiance of the growth chamber. The variance of lamp irradiance was relatively low over time (coefficient of variation, CV ≈ 0.3). Pixel reflectance was calculated by the ratio between pixel radiance and the average of spectralon pixels. We separated reflectance corresponding to vegetation cover pixels from background pixels [58], using a threshold NDVI ≥ 0.6 and we calculated the average of top canopy pixels for each pot (see Figure 2). A radiometric and spectral calibration of the Cubert hyperspectral camera was performed at DTU photonics lab. The spectral response of each pixel was obtained using a monochromator from 430-1030 nm at 1 nm intervals. Subsequently, a radiometric calibration was performed using an integrating sphere to obtain a calibration factor for each pixel in each specific band under nine different illumination levels. The camera integration time was set from 0.5 to 19.9 ms (0.1 ms intervals) at nine different illumination levels. In the growth chamber, we used a spectralon target to characterize the average incoming spectral radiation from the lamps and its fluctuation. On five dates, we acquired a hyperspectral image of the spectralon before measuring the crops. For the remaining dates, we used the average irradiance of the growth chamber. The variance of lamp irradiance was relatively low over time (coefficient of variation, CV ≈ 0.3). Pixel reflectance was calculated by the ratio between pixel radiance and the average of spectralon pixels. We separated reflectance corresponding to vegetation cover pixels from background pixels [58], using a threshold NDVI≥ 0.6 and we calculated the average of top canopy pixels for each pot (see Figure 2).

Thermal Images
The thermal camera was previously calibrated with a Landcal P80P black body radiation source (Land Instruments, Leicester, United Kingdom), presenting an accuracy of 0.95 °C root mean square deviation (RMSD) inside the RERAF at differing ambient and target temperatures. For more details, refer to Köppl et al. 2016 and [65,66]. To extract leaf radiometric temperature (TL,Rad) in each pot, we pre-selected vegetation cover pixels with an automatic procedure, using leaf contour detection for plant recognition, together with histogram thresholding (Figure 2). We then calculated

Thermal Images
The thermal camera was previously calibrated with a Landcal P80P black body radiation source (Land Instruments, Leicester, United Kingdom), presenting an accuracy of 0.95 • C root mean square deviation (RMSD) inside the RERAF at differing ambient and target temperatures. For more details, refer to Köppl et al. 2016 and [65,66]. To extract leaf radiometric temperature (T L,Rad ) in each pot, we pre-selected vegetation cover pixels with an automatic procedure, using leaf contour detection for plant recognition, together with histogram thresholding (Figure 2). We then calculated the average temperature for vegetated pixels. The results of the algorithm were validated with visually photo interpreted leaf temperature from 15 thermal images. The images selected covered a full range of temperature variability. For further details, refer to Gulyas et al. (2020) [67].

Spectral Indices to Predict Leaf Physiology
Due to soil water stress, plants adjust their physiology, morphology, and biochemistry. For instance, on a general basis, under drought, leaves tend to close stomata limiting water vapor loss via transpiration, which reduces photosynthetic carbon uptake and increases leaf temperature. In addition, leaf turgor and chlorophyll content are reduced, affecting growth [13,22]. Note that chlorophyll, carotenoids and anthocyanins provide vital details about plant physiology status, absorbing energy, harvesting light, and protecting leaves from excessive light [68]. Thus, observable changes in leaf temperature and vegetation indices (VI), related to structure, pigments, and photosynthetic efficiency can provide valuable information about how leaf physiology is affected by drought. In our study, we calculated leaf-air temperature difference (∆T) and a total of 13 narrow-band indices ( Table 2), that have previously shown sensitivity to soil water stress. Table 2. Remote sensing indices with equation, definition and reference (ρ = reflectance in specific wavelength, Max dρ = maximum value of the first derivative of reflectance, T L,Rad = leaf radiometric temperature ( • C) and T air = air temperature ( • C)).
Indirect changes on chlorophyll content, related to reduced soil water, can be assessed by calculating the ratio R 700/670 . This VI includes a chlorophyll absorption feature at 670 nm [42,43] and the 700 nm band located at the red-edge region, which is sensitive to chlorophyll concentration [71,75]. VIs tracking chlorophyll content such as the transformed chlorophyll absorption in reflectance index (TCARI) and the ration between TCARI and OSAVI have been successfully used to estimate water stress [10,37]. In addition, TCARI/OSAVI reduces the effect due to changes on LAI and soil reflectance under soil water stress [39,43]. The modified anthocyanin reflectance index (mARI) combine blue, red, and NIR bands to predict changes on anthocyanin content [68]. At 470 nm (blue band), carotenoid absorption can be detected [40] and this wavelength is used in the pigment specific normalized difference index (PSNDc).
The photochemical reflectance index (PRI) is used to estimate photosynthetic activity [41,76]. Because PRI is sensitive to structural changes, pigment levels, soil background, illumination effects, and viewing angles [42,43], other derivations have appeared such as the PRI (570-515) and the sPRI, which have been proved to be more effective to detect water stress [43,76].

Analysis of Variance and Post-Hoc Test
Two-way analysis of variance (ANOVA) was used to investigate the effect of varying soil water and the effect of growth over time on the measured physiological, morphological, and biochemical parameters, and remote sensing derived indices. ANOVA was applied using Python package Pingouin [77]. From this analysis, we obtained the f-value. For better visualization, we also used post-hoc test to perform pair-wise comparisons (t-test) between soil water levels and crops over the seasonal average. Finally, p-value and the coefficient of determination (R 2 ) were obtained to explore the significant relationships between leaf physiology (A, g s , and T r ), SPAD chl and h c with VIs and ∆T. During this study, we applied a significance level of 5%, hence we considered statistical significance at p-value < 0.05. Note that in some cases significance level up to 0.1% was also reported (p-value < 0.001).

Partial Least Squares Regression (PLS-R) Modeling
Partial least squares regression (PLS-R) is recently becoming popular in remote sensing spectroscopy to predict plant traits [15,53,[55][56][57][58][59]61,62,78]. PLS-R is a machine learning technique, typically used to predict an outcome (y) from a set of predictors (x 1 , x 2 , . . . , x p ). In addition, PLS-R can handle highly correlated predictors. In a nutshell, PLS-R decomposes the predictor matrix (X) into a set of loadings and scores with the objective of maximizing co-variance between the scores and y. This is repeated for a given number of latent variables (LV) as the number of loadings and scores necessary to explain sufficient variance in y. The optimal number of LV can be found by employing cross validation [79,80].
We modeled leaf A, leaf g s , and leaf T r using VNIR reflectance, T L,Rad and h c as predictors. In our study, although h c was manually measured, it was defined as a predictor as it can be estimated by photogrammetry techniques with UAVs [52] or from LiDAR sensors. In total, the number of samples totaled 144 for soybean and 126 for maize (Table 1). To assess the synergies of using only hyperspectral or adding other types of remote sensing information, we conducted PLS-R modeling multiple times: (i) with only hyperspectral data (Hyper), (ii) adding T L,Rad , and (iii) including h c , for both crops together, for only soybean and only maize data. PLS-R was implemented using the Python package Scikit-learn [81,82].

PLS-R Preprocessing
Growth and spatial effects are often considered a challenge when using remote sensing and physiology measurements [20]. To remove undesirable obstructive systematic variation, we carried out pre-processing of the spectral data. Standard normal variate (SNV) [83] was applied performing mean and variance normalization for each observation using the equation where N accounts for the number of samples and K for the number of optical hyperspectral wavelengths. µ i and σ i are the mean and standard deviation of spectral data in each sample, respectively. x i,j indicates the value of the i'th number observation in the j'th predictor and x i,j is the corrected predictor after using SNV with zero mean and unit variance. Standardization of the predictor data was subsequently performed (Equation (2)) as suggested by Geladi and Kowalski (1986) [80] and also applied in other remote sensing studies to predict crop function [15,53].
where µ j and σ j are the mean and standard deviation of samples in each spectral wavelength, respectively. = x i,j indicates the corrected predictor after standardization. Similarly, standardization of response data was applied as where µ and σ are the mean and standard deviation of samples for the selected response variable, respectively. y i indicates the corrected response after standardization.

PLS-R Selection and Assessment
The dataset was randomly split into calibration (70%) and validation (30%) sets [55,62]. Using the same dataset, we applied PLS-R over 500 bootstrapped iterations to evaluate the uncertainty of PLS-R models. Bootstrapping is a technique that generates random samples of calibration and validation datasets. Bootstrapping has been found to be the most suitable strategy to validate predictability, ensuring reproducibility and generalizability of the final model [84]. In order to avoid 'over-fitting' and find the optimal number of LV, we applied 10-fold cross validation over the calibration sets by minimizing the corresponding mean squared error over 15 LV. The validation sets were used to estimate generalizability and evaluate the predictive power of the final model. Model performance was assessed using the coefficient of determination (R 2 ), the normalized root mean squared error (NRMSE) (Equation (4)), and the model bias (Equation (5)) of the calibrated and validated datasets, after pooling replicates and the 500 iterations together while separating between crops, dates, and soil water groups.
where y i and y i are the known and predicted values, respectively. S is the total number of samples corresponding to calibration or validation sets. y max and y min represent the maximum and minimum known values. From PLS-R, we assessed the contribution of each predictor to the response using the variable importance of projection (VIP). VIP explains the contribution of individual wavelengths to the overall fitted PLS-R model across components. VIP scores lower than the unit indicate that the contribution of that variable can be neglected. VIP scores were calculated according to Andersen and Bro (2010) [85] as in Equation (6) VIP where J represents the number of variables, F is the total number of components, w j f is the weight of component f and variable j, SSY f is the sum of squares of explained variance for the specific component f and SSY total outlines the total sum of squares explained by the response.
To identify wavelengths correlated to T L,Rad and h c in the latent variable space, we displayed the loadings of the first and second individual variables in that subspace. The goal was to detect which of the 138 VNIR bands provided similar information as T L,Rad and as h c . Figure 3 compares the sensitivity of seasonal canopy evapotranspiration (ET), leaf transpiration rate (T r ), leaf photosynthetic CO 2 assimilation rate (A), and leaf stomatal conductance (g s ) towards changing water availability based on the three water holding capacity (WHC) groups. Comparing among soil water status from 100% (wet/control) to 40%, (dry) soybean showed significant decreases in season average ET, T r , A, and g s (p-value < 0.05). However, for maize the season average responses of leaf A, g s , and T r were not sensitive to soil drought, while they were significantly higher for canopy ET from dry to wet soil water status. Maize A was significantly higher than in soybean plants for the same soil water stress, consistent with higher efficiency of C 4 metabolic path. T r and g s were also higher in maize than soybean, whereas maize ET at canopy level was significantly lower in each water group.

Physiological, Biochemical, and Morphological Responses
Separating the influence of soil drought (WHC) and phenology (days of measurement, DoM) ( Table 3), we saw significant differences for both treatments in soybean physiology (A, T r , and g s ) and biochemistry (chl) at the leaf level, but not for their interaction. Maize leaf parameters only showed a significant response to phenology, with similar variability (F-values) as in soybean parameters.
At canopy level, ET of both crops were sensitive to soil drought and phenology (Table 3 and Figure 3), whereas h c only changed over time (Table 3 and Figure S1). In addition, phenology effects on canopy level parameters (ET and h c ), were more than 150 times higher than replicates variability (F-values > 150) for maize, but not for soybean plants (Table 3). Table 3. Two-way ANOVA test results of season physiological (ET = canopy evapotranspiration, A = photosynthetic CO 2 assimilation rate, T r = transpiration rate, g s = stomatal conductance), biochemical (chl = SPAD chlorophyll content) and morphological parameters (h c = canopy height and LAI = leaf area index) for water holding capacity (WHC) groups and days of measurements (DoM). Numbers indicate F-value and symbols p-value (* p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001 and ns = not significant).  Physiological, biochemical, and morphological parameters ( Figure S3) showed a normal distribution over time among WHC groups and replicates. A strong relation was observed between Tr and gs, showing p-value < 0.001 and coefficients of determination (R 2 ) of 0.75 and 0.86 in soybean and maize, respectively. gs and Tr vs. A for each crop revealed slightly lower relations (R 2 between 0.56 and 0.63, p-value < 0.001). These correlations between physiological traits suggest adequate and regular plant behavior inside the growth chamber ( Figure S3). Coefficients of determination were relatively high between ET and hc for both crops, which is reasonable due to the increases in LAI over time. Maize A, gs, and Tr also correlated with hc, in agreement with the ANOVA (Table 3), showing that variations in leaf physiology were significant across time.

Remote Sensing Responses
For illustrative purposes Figure 4 displays spectra from soybean and maize plants, comparing both crops among soil water levels ( Figure 4A,B) and days of measurements (DoM) ( Figure 4C,D). The seasonal average reflectance of soybean was about 10% higher than maize in near infrared (NIR) bands (760-950 nm), which could indicate higher seasonal LAI for soybean than maize [42]. In the red-edge and NIR regions, reflectance was slightly lower for wet plants (100% WHC) for both crops ( Figure 4A). In the VIS, wet maize plants also showed lower NIR reflectance, while soybean plants Interestingly, comparing wet and dry plants at each day of measurement, soybean and maize ET showed significant decrease since DoM 4. Besides that, soybean leaf chl significantly increased from very early in the growth cycle. However, other parameters such as A, T r , g s or h c were not sensitive in most of the dates ( Figure S2).
Physiological, biochemical, and morphological parameters ( Figure S3) showed a normal distribution over time among WHC groups and replicates. A strong relation was observed between T r and g s , showing p-value < 0.001 and coefficients of determination (R 2 ) of 0.75 and 0.86 in soybean and maize, respectively. g s and T r vs. A for each crop revealed slightly lower relations (R 2 between 0.56 and 0.63, p-value < 0.001). These correlations between physiological traits suggest adequate and regular plant behavior inside the growth chamber ( Figure S3). Coefficients of determination were relatively high between ET and h c for both crops, which is reasonable due to the increases in LAI over time. Maize A, g s , and T r also correlated with h c , in agreement with the ANOVA (Table 3), showing that variations in leaf physiology were significant across time.

Remote Sensing Responses
For illustrative purposes Figure 4 displays spectra from soybean and maize plants, comparing both crops among soil water levels ( Figure 4A,B) and days of measurements (DoM) (Figure 4C,D). The seasonal average reflectance of soybean was about 10% higher than maize in near infrared (NIR) bands (760-950 nm), which could indicate higher seasonal LAI for soybean than maize [42]. In the red-edge and NIR regions, reflectance was slightly lower for wet plants (100% WHC) for both crops ( Figure 4A). In the VIS, wet maize plants also showed lower NIR reflectance, while soybean plants showed the opposite ( Figure 4B). Over time, there was a marked increase in reflectance, especially for maize plants and in NIR wavelengths ( Figure 4C,D).
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 32 showed the opposite ( Figure 4B). Over time, there was a marked increase in reflectance, especially for maize plants and in NIR wavelengths ( Figure 4C,D). Analysis of seasonal responses to soil drought and phenology responses on remote sensing variables is shown in Table 4, Figure 5, and Figure S1. Soybean TL,Rad was significantly lower between dry (40% WHC) and wet (100% WHC) plants, while TL,Rad was not sensitive to water changes in maize ( Figure 5 and Table 4). TL,Rad of both crops was also sensitive to phenology and its interaction with soil water stress (Table 4). Since the air temperature (Tair) was practically constant, ΔT showed similar outcome to TL,Rad for both crops ( Figure S1 and Table 4).
Mean differences due to soil drought for some VIs are displayed in Figure 5, while the rest are in Figure S1. VIs used to predict structural changes such as changes in LAI (e.g., NDVI, EVI, OSAVI, Analysis of seasonal responses to soil drought and phenology responses on remote sensing variables is shown in Table 4, Figure 5, and Figure S1. Soybean T L,Rad was significantly lower between dry (40% WHC) and wet (100% WHC) plants, while T L,Rad was not sensitive to water changes in maize ( Figure 5 and Table 4). T L,Rad of both crops was also sensitive to phenology and its interaction with soil water stress (Table 4). Since the air temperature (T air ) was practically constant, ∆T showed similar outcome to T L,Rad for both crops ( Figure S1 and Table 4). Table 4. Two-way ANOVA test results of seasonal leaf radiometric temperature (T L,Rad ), leaf-air difference temperature (∆T), and vegetation indices (VI), calculated from the hyperspectral sensor described in Table 2, for water holding capacity (WHC) groups and days of measurements (DoM). Numbers indicate F-value and symbols p-value (* p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001, and ns = not significant). Mean differences due to soil drought for some VIs are displayed in Figure 5, while the rest are in Figure S1. VIs used to predict structural changes such as changes in LAI (e.g., NDVI, EVI, OSAVI, and REIP) were statistically significant between wet and dry soil water stratus for maize, but not for soybean. In addition, they changed significantly with phenology for both crops, with higher variation on maize as they showed notably higher F-values compared to other VIs (Table 4). rNDVI, a suitable index used to detect water stress [76], showed significant increase with soil drought ( Figure 5).

Crop
VIs used to detect concentrations of chlorophyll pigments (e.g., TCARI and R (700/670) ), carotenoids (PSNDc), and anthocyanins (mARI), and VIs related to photochemical efficiency (e.g., sPRI and PRI (570-515) ) significantly decreased with soil water stress in soybean plants ( Figure 5 and Figure S1), whereas in maize they showed higher variability along time (Table 4). In addition, it is worth noticing that soybean R (700/670) and sPRI (calculated from 560 and 510 nm) did not significantly change over time (Table 4).
Scatterplots shown in Figure 6 and Figure S4 between VIs with leaf physiology, h c and chl can provide extra information about type of relation, the direction of change, linearity, and saturation of VIs. For instance, both crops presented weak positive linear relations (R 2 < 54%) between VIs with leaf A, T r , and g s (Figure 6 and Figure S4) and negative relationships between ∆T and leaf physiological parameters ( Figure 6).
∆T was not sensitive to changes in maize g s , T r , or A, while it was negatively correlated with soybean g s and T r , with a relatively high correlation for g s (R 2 = 0.47 and p-value < 0.001) ( Figure 6). However, maize g s , T r , or A were sensitive to most VIs, but not for soybean. Highest correlations were shown for maize g s and T r compared to all the structural indices in the following order rNDVI, OSAVI, NDVI, EVI, and REIP. Photosynthetic efficiency indices such as sPRI and PRI (570-515) ; and some VIs associated to pigments like TCARI, PSNDc, and mARI also presented correlations with similar order of magnitude (R 2 between 0.3 and 0.45 with p-value < 0.01) (Figure 6 and Figure S4).    Table 2) and leaf-air temperature difference (ΔT in °C) vs. leaf photosynthetic CO2 assimilation rate (A in µmol CO2 m −2 s -1 ), leaf stomatal conductance (gs in mol H2O m −2 s −1 ), leaf transpiration rate (Tr mmol H2O m −2 s −1 ), SPAD chlorophyll content (chl in µg m -2 ), and canopy height (hc in m) for soybean (blue crosses) and maize (red dots). Points represent the average of six replicas. Coefficient of determination (R 2 ) is shown and p-value is represented with stars (* p-value < 0.05; ** p-value <0.01; *** p-value < 0.001).  Table 2) and leaf-air temperature difference (∆T in • C) vs. leaf photosynthetic CO 2 assimilation rate (A in µmol CO 2 m −2 s −1 ), leaf stomatal conductance (g s in mol H 2 O m −2 s −1 ), leaf transpiration rate (T r mmol H 2 O m −2 s −1 ), SPAD chlorophyll content (chl in µg m −2 ), and canopy height (h c in m) for soybean (blue crosses) and maize (red dots). Points represent the average of six replicas. Coefficient of determination (R 2 ) is shown and p-value is represented with stars (* p-value < 0.05; ** p-value < 0.01; *** p-value < 0.001).
Soybean h c was positively correlated with most of the VIs related with vegetation structure such as REIP, EVI, and OSAVI (R 2 between 0.4 and 0.6) and some chlorophyll indices like TCARI and TCARI/OSAVI (R 2~0 .35, p-value < 0.01), while negatively correlated to PRI changes (R 2 = 0.54). Maize h c presented higher sensitivity than soybean h c to VIs with most R 2 between 0.73 and 0.96 ( Figure 6 and Figure S4). Note that TCARI and TCARI/OSAVI decrease are linked to chl increase [39,43].
In soybean, the strongest relations with VIs were found for chl with VIs used to detect pigment and photosynthetic efficiency, while no correlations were found with any of the structural VIs, except for rNDVI, which is also associated to water stress and REIP, which can estimate chlorophyll concentration. Even though maize VIs were also sensitivity to chl, the relation was opposite from soybean (e.g., TCARI), same as for R 700/670 or PRI. Maize VIs related to structural changes were sensitive to most of the variables (chl, h c , T r , g s , and A), which was not the case for soybean (e.g., NDVI, EVI) ( Figure 6 and Figure S4).

Partial Least Squares Regression (PLS-R) to Predict Leaf Physiology
Before using PLS-R, a SNV pre-processing was applied to minimize growth and crops differences ( In soybean, the strongest relations with VIs were found for chl with VIs used to detect pigment and photosynthetic efficiency, while no correlations were found with any of the structural VIs, except for rNDVI, which is also associated to water stress and REIP, which can estimate chlorophyll concentration. Even though maize VIs were also sensitivity to chl, the relation was opposite from soybean (e.g., TCARI), same as for R700/670 or PRI. Maize VIs related to structural changes were sensitive to most of the variables (chl, hc, Tr, gs, and A), which was not the case for soybean (e.g., NDVI, EVI) ( Figure 6 and Figure S4).

Partial Least Squares Regression (PLS-R) to Predict Leaf Physiology
Before using PLS-R, a SNV pre-processing was applied to minimize growth and crops differences (Figure 7 compared to Figure 4).  Table 5 provides R 2 , NRMSE and bias from PLS-R to assess model performance differences when using only hyperspectral data or when including TL,Rad and hc, pooling data from both crops together or separated. Calibration NRMSE of 14% and below as well as test set NRMSE of 18% and below indicated good model performance with 500 permutations. In addition, all the R 2 values were statistically significant (p < 0.001). Table 5. Fit statistics for 500 bootstrapped PLS-R model for each crop and both crops (soybean + maize) together when applying PLS-R for hyperspectral data (Hyper) only, including leaf radiometric temperature (TL,Rad), and finally adding canopy height (hc). The coefficient of determination (R 2 ), normalized root mean squared error (NRMSE) and bias are provided for PLS-R Calibration fitting (cal) during the PLS-R model development with 70% of the data and for PLS-R validation (val) with 30% of the data. These metrics are calculated from bootstrapped data for the same crop, water group, date and set. LV indicates the average of latent variables over 500 iterations, after using 10-fold cross validation for each iteration.   Table 5 provides R 2 , NRMSE and bias from PLS-R to assess model performance differences when using only hyperspectral data or when including T L,Rad and h c , pooling data from both crops together or separated. Calibration NRMSE of 14% and below as well as test set NRMSE of 18% and below indicated good model performance with 500 permutations. In addition, all the R 2 values were statistically significant (p < 0.001).
It is notable that the best model performance (both R 2 and NRMSE) for the three physiological variables studied (A, g s , and T r ) was found for joint crop models rather than separate crop models ( Table 5). This joint crop model (Figure 8), showed higher correlations between predicted and observed data, and lower calibration and validation errors, especially for A (R val 2 = 0.92 and NRMSE val = 8%), followed by T r (R val 2 = 0.87 and NRMSE val = 11%), and g s (R val 2 = 0.82, NRMSE val = 11%). This is most likely due to the larger range of the leaf physiological variables and the increased number of observations ( Figure S3 and Figure 8). As previously shown in Figure 3, maize A, g s , and T r presented higher response values than soybean. Matthes et al. (2015) [55] also observed better PLS-R scores in pasture and rise gross primary productivity (GPP) and net ecosystem exchange (NEE) over Mediterranean climate, using hyperspectral data (400-900 nm). Comparing between crops, maize presented better predictions for T r and A while g s was better predicted in soybean ( Table 5). The main improvement from adding T L,Rad to PLS-R models was found in soybean models (Table 5), especially for g s . Soybean g s model including T L,Rad , increased about 25% R 2 val (13% R 2 cal ), decreased 13% NRMSE val (16% NRMSE cal ), and reduced 10% bias val (13% bias cal ) compared to just using hyperspectral data. Adding h c was mostly improving maize models, especially T r , reducing NRMSE val about 15% and bias val 30%. Maize g s and T r models showed increased R 2 calibration and validation values of about 10%.
VIP scores (Figure 9) higher than 1 indicate important waveband regions to predict A, g s , and T r . Overall, we observed the same bands contribution for soybean A, g s , and T r , while in maize, the role of specific narrow bands, was not so clear, presenting lower VIP scores and similar influence over wavelengths. In maize A, gs, and Tr, we could highlight some important regions around the red edge and NIR, more specifically, at the beginning and end of the red edge slope with peaks around 700 and 770 nm. In the NIR of maize A VIP, some displaced regions were found compared to gs and Tr, which showed a peak at 843 nm. In addition, further in the NIR, different wavelengths were prominent for physiology, with some relations between A with gs and Tr. Maize gs and Tr VIP scores were also influenced by green bands (520 and 560 nm).
For A, gs, and Tr in the studied crops, unimportant regions were observed in blue and yellowred bands from 450 to 500 nm and between 580 and 680 nm, respectively. In addition, VIPs for soybean and maize gs and Tr were very low around 900 nm. Based on VIP threshold score and loading weights (Figure 9), the most influential regions for soybean A, g s , and T r were the green band in the VIS (peak around 516 nm), at the beginning (691 nm) and end (780-783 nm) of the red-edge slope and in the NIR plateau region (peaks at 831 and 850 nm). Furthermore, g s and T r showed a spike in another green band (540 nm), especially important in g s VIP curve, presenting the second highest value. In soybean g s and T r , the highest VIP score was given by T L,Rad contribution. In maize, however, h c played more essential role as shown in Table 5. For A of both crops, these two added predictors did not seem to add much information. Note that h c was also important for predicting soybean T r . To assess the correlation between hyperspectral, thermal, and hc data, we plotted the predictors PLS-R mean loadings of the 500 bootstrapped iterations for the first and second LVs (Figure 10), highlighting important wavelengths from the VIP analysis ( Figure 9). From Figure 10, we observed that most of the noteworthy bands from VIP analysis were positively or negatively correlated with hc and TL,Rad for both crops, especially for soybean plants.
From Figure 10, In soybean A, gs, and Tr, we found that wavelengths corresponding to absorption of chlorophyll b (516 nm) and NIR plateau (850 nm) were the closest to TL,Rad and in the same quadrant. Opposite trends to TL,Rad were shown by red-edge (780-783 nm) and NIR canopy moisture (831 nm) bands. Similarly occurred with the green band at 540 nm for soybean gs and Tr. Interestingly in three soybean physiological parameters, the chlorophyll a absorption band (691 nm) was near to hc.
In maize A, gs, and Tr, wavelengths at the end of the red edge slope (760-780 nm) were the most related to hc. Also correlated with hc but with negative LV1, bands at the beginning of the red-edge (698-706 nm) and NIR regions (930-950 nm were highlighted. Note that in accordance with previous findings, TL,Rad was not considerably contributing to maize physiological modeling. In maize A, g s , and T r , we could highlight some important regions around the red edge and NIR, more specifically, at the beginning and end of the red edge slope with peaks around 700 and 770 nm. In the NIR of maize A VIP, some displaced regions were found compared to g s and T r , which showed a peak at 843 nm. In addition, further in the NIR, different wavelengths were prominent for physiology, with some relations between A with g s and T r . Maize g s and T r VIP scores were also influenced by green bands (520 and 560 nm).
For A, g s , and T r in the studied crops, unimportant regions were observed in blue and yellow-red bands from 450 to 500 nm and between 580 and 680 nm, respectively. In addition, VIPs for soybean and maize g s and T r were very low around 900 nm.
To assess the correlation between hyperspectral, thermal, and h c data, we plotted the predictors PLS-R mean loadings of the 500 bootstrapped iterations for the first and second LVs (Figure 10), highlighting important wavelengths from the VIP analysis ( Figure 9). From Figure 10, we observed that most of the noteworthy bands from VIP analysis were positively or negatively correlated with h c and T L,Rad for both crops, especially for soybean plants.
From Figure 10, In soybean A, g s , and T r , we found that wavelengths corresponding to absorption of chlorophyll b (516 nm) and NIR plateau (850 nm) were the closest to T L,Rad and in the same quadrant. Opposite trends to T L,Rad were shown by red-edge (780-783 nm) and NIR canopy moisture (831 nm) bands. Similarly occurred with the green band at 540 nm for soybean g s and T r . Interestingly in three soybean physiological parameters, the chlorophyll a absorption band (691 nm) was near to h c .
In maize A, g s , and T r , wavelengths at the end of the red edge slope (760-780 nm) were the most related to h c . Also correlated with h c but with negative LV1, bands at the beginning of the red-edge (698-706 nm) and NIR regions (930-950 nm were highlighted. Note that in accordance with previous findings, T L,Rad was not considerably contributing to maize physiological modeling. Figure 10. PLS-R mean loadings of the 500 bootstrapped iterations for the predictor variables: VNIR wavelengths (color bar from blue to red colors with increasing nm), TL,Rad (black star), and hc (black triangle) for the first latent variable (LV, x-axis) and the second LV (y-axis) comparing photosynthetic CO2 assimilation rate (A), leaf stomatal conductance (gs), and transpiration rate (Tr) of soybean (left figures) and maize (right figures). Important wavelengths from VIP analysis are marked with a black dot and the corresponding nm number. Wavelengths near to TL,Rad and hc as well as wavelengths in the opposite quadrant to these two predictors have similar explanator power. Figure 10. PLS-R mean loadings of the 500 bootstrapped iterations for the predictor variables: VNIR wavelengths (color bar from blue to red colors with increasing nm), T L,Rad (black star), and h c (black triangle) for the first latent variable (LV, x-axis) and the second LV (y-axis) comparing photosynthetic CO 2 assimilation rate (A), leaf stomatal conductance (g s ), and transpiration rate (T r ) of soybean (left figures) and maize (right figures). Important wavelengths from VIP analysis are marked with a black dot and the corresponding nm number. Wavelengths near to T L,Rad and h c as well as wavelengths in the opposite quadrant to these two predictors have similar explanator power.

Physiological, Biochemical, Morphological, and Remote Sensing Responses under Water Stress
Changes in leaf physiology are the result of plant-level regulation, and therefore changes in g s , T r , and A associated with soil drought have to be interpreted by looking at transpiration responses of the overall canopy. After 4-7 DoM, leaves were practically covering the full soil surface, especially in soybean plants. Thus, soil evaporation could be neglected, assuming canopy transpiration similar to ET. Our data showed that for a 60% soil water reduction, both crops experienced a significant reduction about the same order of magnitude in average canopy transpiration (around 23% for maize and 35% for soybean) (Figure 3). However, maize and soybean plants presented different responses in T r and g s at the leaf level.
For soybean, a proportional decline of leaf g s (≈34%) and leaf T r (≈24%) was observed after the reduction of 35% ET (Figure 3). This suggests that the strategy of soybean plants to cope with drought is through tight stomatal control, also shown by the significant rise of leaf vapor pressure deficit (VPD leaf ≈ 13%) ( Figure S1) and T L,Rad (≈3%) (Figure 3). From this study, we found that T L,Rad is a key remote sensing index to track drought responses in soybean plants, as it captures the warming effect from less cooling from lower T r and g s (Figure 3). Furthermore, ∆T was the only index significantly correlated with soybean g s and T r ( Figure 6). Similar response in g s of soybean plants has been also shown by others [17,19,86]. Thigh stomatal control could be associated to soybean plants not acclimating well to soil water stress due to lack of adjustment in leaf hydraulic conductance [87]. Thus, to avoid hydraulic failure, it tends to close stomata. As in soybean plants, most of the reduction in ET can be linked to reductions in leaf T r , water stress impacts on growth (i.e., LAI decreases) did not seem to play a big role ( Table 3). The main support for this argument relies on LAI measurements at the end of the season with not significant differences between wet and dry soil water levels ( Table 3 and Figure S5). In addition, remote sensing indices traditionally associated with LAI such as NDVI [42,53,88], OSAVI [37,39,85], and REIP [75,88], were not significantly affected by changes in soil water (Table 4, Figure 5 and Figure S1). It is possible that, as VIs were calculated only from vegetation pixels (Figure 2), their ability to detect decrease in LAI might be reduced. However, it can be seen that h c , a proxy for cell enlargement and associated to growth [7,17], did not suffer significant changes due to drought (Table 3). Because VIs were obtained from average pixels of top vegetation cover from a high resolution camera (10.5 mm/px), we assumed that most of the detected reflectance by the sensor was coming from top of canopy leaves and it is not considered as a composite plant-canopy reflectance. However, some of the pixels might incorporate reflectance of several overlapping leaves ( Figure 2) and leaf scattering effects might be observed on the NIR part of the spectrum [32], which can be also related with changes in leaf angle due to drought effects [10,32,89]. By removing soil background pixels (Figure 2), we eliminated most of the reflectivity from soil, supported by similar behavior between NDVI and OSAVI, which minimizes soil effects [70], with significant differences for soybean but not for maize (Table 3). Thus, we could infer that soil brightness effects were low. Nevertheless, for further improvements on reflectance measurements, 'bare soil' reflectance for each soil water level can provide useful information to reduce soil brightness effects [90].
In maize, all the leaf physiological variables (g s , T r , and A) (Figure 3) as well as VPD leaf and intercellular CO 2 concentration (C i ) ( Figure S1) were very similar among WHC groups. Nonetheless, such differences in T r and A were enough to result in a significantly lower water use efficiency (WUE ≈ 10%) with 60% reduced soil water availability ( Figure S1). Therefore, the detected 23% drop of ET under soil water stress (Figure 3) must be mostly driven by a reduction in LAI, as morphological adjustments are an effective mechanism to cope with water stress in maize plants [8].
Our results support this hypothesis, where measurements of LAI performed at the end of the experiment were significantly lower with water stress for maize but not for soybean (Table 3 and Figure S5). Moreover, structural related VIs such as NDVI, rNDVI, OSAVI, EVI, and REIP showed significant effects between soil water levels ( Table 4) and significantly positive relations between these VIs with g s and T r (Figure 6 and Figure S4). In recent studies, structural VIs have also successfully been used to predicted water stress proxies (e.g., g s , leaf water content, and leaf water potential) in maize plants [91,92]. In our study, the effect of growth through these VIs is shown by the high sensitivity of phenology (DoM), lack of significance with soil drought (WHC) ( Table 4), and the strong positive correlations with h c (R 2 > 0.86) (Figure 6 and Figure S4). In addition, maize reflectance presented higher reflectance over time than soybean ( Figure 4C,D), which indeed could imply more variability due to growth, as changes in LAI can be detected in the NIR region close to the red-edge [56,75]. In contrast to soybean, maize T L,Rad was not sensitive to drought (Table 4) and ∆T did not vary with phenology. In addition, ∆T was the only remote sensing variable not related to any maize parameter (g s , T r , A, h c , and chl) ( Figure 6). Unlike our results, others studies of greenhouse maize plants have been successful in detecting drought through thermal imaging, likely due to higher air temperature (27-28 • C) increasing soil water stress [50,93].
Despite of such a large mean reduction in soybean g s from wet (100%) to dry (40%) WHC, soybean plants significantly increased about 15% their WUE ( Figure S1) as they only experienced a small decrease in A (≈12%) relative to the reduction in transpiration at canopy and leaf levels ( Figure 2). These could be explained by a 16% mean increase in top leaves chl concentration, measured with SPAD ( Figure S1). Zhang et al. (2016) [19] also found increased chlorophyll content as well as decreased g s , T r , and A in greenhouse soybean plants with soil water reductions of about 40%. However, they reported declines in WUE of about 45%. Wijewardana et al. (2019) [17] have interpreted WUE increases as an adaptive mechanism to drought due to transpiration decline under soil moisture stress in soybean plants. However, in their experiment, they observed a 24% decrease on chlorophyll content of soybean plants with 38% increase in carotenoids. Explanations to differences among studies might be associated to experimental design. For instance, plants in Wijewardana et al. (2019) [17] were exposed to higher radiation and temperatures (29 • C) than in our experiment, with a longer cycle, which would increase soil water stress. Cultivar types also present different morphological and physiological adjustments to drought [7,86].
Moreover, VIs used to detect changes in chlorophyll concentration (TCARI, TCARI/OSAVI, and R 700/670 ) can track responses to soil drought across different dates, especially for R 700/670 (Table 4, Figure 6 and Figure S4). Unlike soybean, maize leaf chl did not change across soil water levels (Table 3), reflected also in the lack of significance in TCARI and TCARI/OSAVI indices. However, this was not the case for R 700/670 (p-value < 0.01) ( Table 4), likely because this VI does not minimize LAI variations as TCARI and TCARI/OSAVI [39]. Differences on chl response between crops can be related to decreased soybean reflectance in the VIS with soil drought (up to around 640 nm), whereas maize reflectance increased ( Figure 4A,B). Similarly, Feng et al. (2013) [94] found greater maize reflectance under soil water stress in this region. VIs linked to photosynthetic efficiency like PRI are often used to detect soil moisture changes [37,41,43,45,76]. However, PRI is also quite sensitive to structural changes, pigment levels, soil background, illumination effects, and viewing angles [42,43]. As the experiment was conducted in a growth chamber with high-pressure mercury and halogen lamps (peak at 538 nm), some of the PRI results might not be like in field conditions ( Figure 3). However, variations of PRI calculated with other bands as sPRI (510 nm and 560 nm) and PRI 570-515 were significantly reduced presenting less photosynthetic activity with drought (Table 4), especially for soybean plants (C 3 ) as they are less efficient than C 4 species (maize) [21,22]. Note that soybean sPRI as well as R 700/670 were not influenced by changes across time (Table 4), making them attractive to detect drought in soybean plants, likely subtracting growing effects. Unlike soybean, not a single maize VI was only significant across water groups (Table 4).
One area of interest is the possibility to use remote sensing to detect isohydric and anisohydric behavior, which is concept related to leaf water potential (Ψ L ) [95]. Leaf water content can be retrieved from optical (between 800 nm and 2500 nm) and microwave radiation [10]. In our experiment, because our hyperspectral camera only covers the VNIR region (400-900 nm), we did not measure Ψ L . However, if reductions in g s are linked to reductions in Ψ L as shown by Wijewardana et al. (2019) [17], it could be inferred that soybean plants behaved more as isohydric crop compared to our type of maize plants. This could also be associated with the type of photosynthetic pathway as C 3 species such as soybean tend to exhibit more considerable reductions in g s than their C 4 relatives under drought [23]. Nevertheless, according to Costa et al. (2013) [48], g s is better indicator of soil water stress in isohydric plants, enabling thermal imaging to better detect drought, as we observed in our results ( Figure 6). In addition, Konings et al. (2017) [25] also identified a maize cultivar as anisohydric. Still, the concept of iso/anisohydric appears more complex than can be represented as a binary factor, and could instead occur as a response to environmental thresholds or occur along a continuum [95,96]. Therefore, discrepancies in this classification are found for soybean [48,97] and maize crops [16,48,[97][98][99].

Synergies of Optical, Thermal, and Canopy Height Observations for PLS-R Modeling
In response to drought severity and duration, crops undertake different adjustments in plant physiology, morphology, and biochemistry. Some of those changes are noticeable from the drought onset (fast changes), while others are noticeable only after some time, also referred as slow changes [10]. For example, stomatal closure occurs practically instantaneous under water stress, while effects on chlorophyll content and structure tend to be visible after longer periods of time [10,11,13]. Our study showed that for soybean plants T L,Rad was a faster indicator of physiological adjustments than VIS indices such as NDVI and PRI ( Figure 6 and Figure S4) as shown also by [45,100,101]). Nonetheless, improved detection of water stress can be achieved by combining thermal and optical imaging [37,46,47,102], as the energy balance is affected by drought. The temperature is the result of the balancing of several processes, including transpiration cooling and absorption of radiation. For example, maize plants were warmer than soybean plants ( Figure 5), despite presenting higher T r for any soil moisture condition (Figure 3). Figure 4A shows how maize wet plants reflected less energy than wet soybean in the VNIR wavelengths, indicating more energy absorption, translated in significantly higher thermal emission. From the point of view of plant physiology, although h c is a parameter at the canopy level, it integrates leaf level processes in response to stress: e.g., cell enlargement has been shown very sensitive to water stress [11,13], and h c reduction could be associated with decrease on cell enlargement [7,17]. For this reason, h c was used as predictor for water stress. In addition, h c is a growth indicator easily obtained from remote sensing [52], that provided valuable information to model g s , T r , and A ( Figure 9). PLS-R is not only a tool to model crop physiology parameters, but the VIP scores can inform of the more relevant wavelengths and variables responding to soil drought. In fact, the highest VIP scores for g s and T r correspond to T L,Rad in soybean and h c in maize (Figure 9), reflecting the dominant effects of soil drought: stomatal closure on soybean and plant growth on maize plants. Moreover, model improvements were appreciated in soybean g s model when including T L,Rad predictor, and in maize T r model when adding h c (Table 5). A curious fact when looking at hyperspectral band contribution was that soybean crops shared importance on similar wavelengths among A, g s , and T r , whereas maize response varied between traits with low VIP signal (Figure 8). From this, we can speculate that soybean processes are more interdependent than maize mechanisms.
Because part of the information obtained from the new variables (T L,Rad and h c ) could be implicit in some of the 138 bands of the hyperspectral sensor, we looked into the relation between LV1 and LV2 for all predictors ( Figure 10). Overall, we found that most of the bands showing high VIP scores ( Figure 9) were correlated with T L,Rad and h c (Figure 10). These correlations could explain why model performance did not improve remarkably by adding the two new predictors (Table 5). For soybean g s , T r , and A peaks with VIP scores above one were found in regions that correspond to absorption of chlorophyll b (516 nm), chlorophyll a (691 nm), red edge (780-783 nm), and the NIR plateau (831 and 850 nm) (Figure 9). Most of these bands clearly overlap the ones used in the VIs studied ( Table 2), but they are not exactly centered at the same wavelength, highlighting the importance of using hyperspectral remote sensing with contiguous narrowband. Further research would be needed to confirm this finding across different illumination sources and conditions.
Regions around 515 nm and 700 nm have been reported to be very suitable to estimate A and indirectly g s due to their ability to detect chlorophyll changes [42,75]. Similar to our result for soybean, Matthes et al. (2015) [55] found maximum VIP scores between 670 and 680 nm when using PLS-R to predict GPP and NEE in pasture and rise sites under Mediterranean climate. In addition, a recent study concluded, based on PLS-R VIP analysis, that most important bands to predict wheat g s , T r , and A from hyperspectral data were mainly placed in the VIS and red-edge regions [57]. Wavelengths around 780 nm presented one of the highest VIP values in soybean g s , T r , and A, in agreement with mARI being sensitive to T r . Note that narrow wavebands (831 and 850 nm), often combined with water absorption bands in SWIR to obtain water and moisture indices [42], were not used in VI calculations. In Figure 9, soybean T r and g s VIP showed also a spike centered at the green band 540 nm, which might be partly due to artifacts caused by artificial illumination sources ( Figure 4). Nonetheless, Rapaport et al. (2015) [15] also found that this band was one of the bands contributing the most to estimate leaf water potential to detect water stress under solar irradiance conditions.
In relation to artifacts caused by the lamps (Figure 4), to some extent, they should be removed when calculating reflectance by the normalization with incoming irradiance and the normalization prior to PLSR processing. Even though it is possible that some of these artifacts could be related to temporal changes in the lamp's irradiance caused by voltage oscillations, other authors have found similar effects in reflectance due to lamps irradiance spikes. Schuerger & Richards (2006) [103] investigated reflectance of healthy and water stressed pepper plants under seven artificial illumination sources and they found that high pressure sodium and metal halide lamps, with irradiance spikes similar to those in our high pressure mercury and halogen lamps, introduced spikes in the reflectance spectra. Thus, peaks in the spectral irradiance distribution of the lamps, compared to the solar irradiance, may modify mostly plants reflection while the absorbed light might not be affected in the same proportion, and, therefore, plants release the extra light. Nevertheless, these small reflectance artifacts do not seem to show an inconvenient to effectively detect plant changes in reflectance due to water stress under these type of lamps [103].
Unlike soybean, the contribution of wavelengths to maize g s , T r , and A PLS-R modeling was more spread over the whole spectra with lower VIP scores ( Figure 9). Nonetheless, there were a few common regions with spikes at two red-edge wavelengths (around 700 and 760-780 nm). Note that bands around 770 nm were near to h c , explaining similar trend. This was not the case for 700 nm presenting negative loading in the first LV ( Figure 10). These bands are closely related to the rNDVI (705 and 750 nm), which relate to structural changes to detect drought [76]. Bands at 705 nm and 700-740 nm relate to water stress while 760 nm accounts for water absorption [42]. In line with this and using PLS-R, strong relations were found between red-edge bands and V cmax (correlated to A) in nine agricultural crops [58]. In Figure 9, green bands (around 520 and 560 nm) for predicting maize g s and T r showed important contribution, but not same information as h c (Figure 10). These bands were not used in VI calculations, despite being close to PRI and its variants (Table 2). Thus, they could be potential wavelengths to predict maize g s and T r under water stress. NIR band (843 nm) was also important (Figure 9). h c has been previously related with near infrared bands [42], associated with LAI and leaf scattering [32]. Maize A VIP was substantially higher at around 890 and 920-950 nm, while T r VIP was prominent at 952 nm and g s VIP at 935 nm ( Figure 9). In the NIR spectral bands, water absorption features are identified around 900 and 970 nm (water index) [104]. Although our hyperspectral sensor detected changes related to water absorption between 870 nm and 960 nm, these wavelengths were more noisy than VIS bands with wider full width at half maximum ( Figure 9). Thus, the models might not detect high signals in water bands. According to Taylor et al. (2011) [23], in C 4 species, drought effects in A and T r might be more visible in the whole plant rather in specific leaf compared with C 3 plants.
Finally, it is worth mentioning that this study was conducted during a short vegetative period of two crops under specific controlled environmental conditions, presenting certain limitations. Thus, to assure the reproducibility of these results, further research will be required in more species, repeating and testing the method in a larger dataset as recommended by Matthes et al. (2015) [55]. In addition, a step further will be to apply this study under external ambient conditions as similarly done by Rapaport et al. (2015) [15].

Conclusions
In this study, we showed that changes in the function of soybean and maize under water stress are expressed in various reflectance and emission traits, according to their different photosynthetic pathways (C 3 and C 4 ) and water use strategies. In a controlled environment, we demonstrated the joint power of optical hyperspectral (VNIR), thermal (T L,Rad ), and canopy height (h c ), to predict three key physiological variables regulating water and carbon assimilation: leaf stomatal conductance (g s ), leaf transpiration (T r ), and leaf photosynthesis (A) during a vegetative growth season.
Under a 60% reduction in soil water content, canopy transpiration of both soybean and maize decreased significantly (35% and 23%, respectively). In soybean (C 3 ), the reduction of canopy transpiration was achieved by a tighter leaf stomatal control, suggesting a more isohydric behavior, while for maize (C 4 ), gas exchange parameters were not significantly affected, which could indicate a more anisohydric tendency. For maize, the reduction of canopy transpiration was caused by structural changes, with leaf area index (LAI) and h c significantly reduced by the end of the season. These very different responses were readily picked by some narrowband vegetation indices (VIs) and the leaf-air temperature difference (∆T), but only after drought responses were separated from variations in phenology.
Although well-known VIs are useful to point out some of the plants potential mechanisms to cope with water stress, almost none of the studied VIs just responded to water stress. In addition, the VI variance explained individually in relation to water stress was not enough to estimate leaf gas exchange processes. Machine learning methods such PLS-R enable identification and exploitation of other wavelengths and variables without any prior assumption about the exact mechanisms involved in order to predict the studied physiological parameters. However, PLS-R requires a large number of samples and range of variation, which could explain why a joint crop model (R 2 val ≈ 0.77-0.92, NRMSE val < 13%) performed better than individual crop models (R 2 val ≈ 0.54-0.81, NRMSE val < 18%), even though maize and soybean presented quite distinct functioning for the sole purpose of predicting leaf gas exchange. Nonetheless, separated PLS-R crop models were useful to identify the most relevant variables in each crop, such as T L,Rad for soybean and h c for maize according to our experimental conditions. Interestingly, the most important spectral bands sensitive to drought, derived from PLS-R analysis, were not exactly centered at the same wavelengths of the studied VIs sensitive to drought, highlighting the benefit of having contiguous narrow spectral bands to predict leaf physiology and suggesting different wavelengths combinations based on crop type. Additionally, PLS-R VIP results showed that soybean relevant wavelengths were similar for A, g s , and T r , whereas maize response varied between physiological parameters with lower VIP scores, suggesting higher interdependence of processes in soybean than maize.
The results from this study are only a first but a promising step towards larger scale applications, using airborne or satellite remote sensing. PLS-R estimates of leaf-level parameters such as g s could help to parameterize canopy level GPP and ET models, recognize different photosynthetic paths, or identify the degree of stomatal closure in response to drought, which could be linked to isohydric/anisohydric behavior.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/12/19/3182/s1, Figure S1. Comparison between soil water levels: Wet/control (100% WHC), Mid (70% WHC) and Dry (40% WHC) for each crop: soybean (blue boxes) and maize (dashed reddish boxes) h c , chl, C i , WUE, leaf temperature from LI-6400 (T leaf ), VPD Lea f, ∆T, radiometric leaf standard deviation temperature (T L,Rad,std ), OSAVI, mARI, TCARI/OSAVI, PSNDc, PRI (570-515) , and PRI. Dots in maize and crosses in soybean represent the seasonal average of the six replicates for each crop. Boxplots contain six replicates. Boxes show the 25th and 75th percentiles of the interquartile range (IQR), horizontal line represents the median of the data and whiskers extend from the edges of box to show 1.96 IQR. Points outside this range are outliers. Different letters indicate that differences are significant at p-value < 0.05 level. Figure S2. Time series plot per day of measurement (DoM) of wet (100% WHC) and dry (40% WHC) soybean (bluish lines in left figures) and maize (reddish lines in right figures) ET, T r , g s , A, h c , C i , chl, and VPD leaf . Asterisks (*) indicate significant differences between wet (control) and dry plants at 5% significance level. Figure S3. Scatter plots between different variables including physiology (canopy evapotranspiration = ET in mm day −1 , leaf photosynthesis = A in µmol CO 2 m −2 s −1 , leaf stomatal conductance = g s in mol H 2 O m −2 s −1 and leaf transpiration = T r in mmol H 2 O m −2 s −1 ), biochemistry at leaf level (chlorophyll content = chl in µg m −2 ) and morphology (canopy height = h c in m). Soybean (blue crosses) and maize (orange dots) are separated with their respective coefficient of determination (R 2 ) and p-value indicated with stars (* p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001). The diagonal shows the distribution of each parameter for both crops. Figure S4. Linear correlation between vegetation indices (VIs) ( Table 2) vs. A (µmol CO 2 m −2 s −1 ), g s (mol H 2 O m −2 s −1 ), Tr (mmol H 2 O m −2 s −1 ) and chl (µg m −2 ) at leaf level and h c (m) for soybean (blue crosses) and maize (red dots). Points represent the average of six replicas. Coefficient of determination (R 2 ) is shown and p-value is represented with stars (* p-value < 0.05; ** p-value < 0.01; *** p-value < 0.001). Figure   Acknowledgments: The authors would like to thank Pablo Miguel Martinez Sanchez, Liang Shen, and Gorka Mendiguren Gonzalez for the help on experimental set up and measurements, and Hongxiao Jin for the support regarding camera calibrations. The authors would also like to thank Maica Martinez and to Davide Ferrari (SAGEA Centro di Saggio s.r.l, Italy) for the support with the seeds and to DTU RERAF facilities for the support with the growth chamber. In addition, the authors would like to thank the editors and the three anonymous reviewers for the comments and suggestions to improve the quality of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.