Estimating Vertical Distribution of Leaf Water Content within Wheat Canopies after Head Emergence

: Monitoring vertical proﬁle of leaf water content (LWC) within wheat canopies after head emergence is vital signiﬁcant for increasing crop yield. However, the estimation of vertical distribution of LWC from remote sensing data is still challenging due to the effects of wheat spikes and the efﬁcacy of sensor measurement from the nadir direction. Using two-year ﬁeld experiments with different growth stages after head emergence, N rates, wheat cultivars, we investigated the vertical distribution of LWC within canopies, the changes of canopy reﬂectance after spikes removal, the relationship between spectral indices and LWC in the upper-, middle- and bottom-layer. The interrelationship among vertical LWC were constructed, and four ratio of reﬂectance difference (RRD) type of indices were proposed based on the published WI and NDWSI indices to determine vertical distribution of LWC. The results indicated a bell shape distribution of LWC in wheat plants with the highest value appeared at the middle layer, and signiﬁcant linear correlations between middle-LWC vs. upper-LWC and middle-LWC vs. bottom-LWC (r ≥ 0.92) were identiﬁed. The effects of wheat spikes on spectral reﬂectance mainly occurred in near infrared to shortwave infrared regions, which then decreased the accuracy of LWC estimation. Spectral indices at the middle layer outperformed the other two layers in LWC assessment and were less susceptible to wheat spikes effects, in particular, the newly proposed narrow-band WI-4 and NDWSI-4 indices exhibited great potential in tracking the changes of middle-LWC (R 2 = 0.82 and 0.84, respectively). By taking into account the effects of wheat spikes and the interrelationship of vertical LWC within canopies, an indirect induction strategy was developed for modeling the upper-LWC and bottom-LWC. It was found that the indirect induction models based on the WI-4 and NDWSI-4 indices were more effective than the models obtained from conventional direct estimation method, with R 2 of 0.78 and 0.81 for the upper-LWC estimation, and 0.75 and 0.74 for the bottom-LWC estimation, respectively.


Introduction
Leaf water content (LWC) in crops is closely associated with physiological processes and morphological structure of plants [1]. In agricultural practice, proper application of fertilizer can promote the growth of plant root, then enhance the ability of water uptakes from soil [2,3]. Abiotic stresses (e.g., nitrogen or phosphorus stress), by contrast, generally reduce crop water conductivity and lead to a decreased LWC [4], which would slow the rate of photosynthetic carbon assimilation and restrain the productivity and growth of crops [5,6]. Therefore, accurate estimation and real time detection of LWC and its vertical distribution with crop canopies, particularly in various functional leaf layers, is helpful for understanding the real water status and water transfer process in plants, and making sound decisions for effective water and fertilization management in precision agriculture.
Remote sensing technology has been successfully applied in monitoring physical and chemical changes in crops across different scales [7][8][9][10]. The fundamental of remote estimation is the interaction between plant chemical compositions, such as LWC, with electromagnetic radiation [11]. Numerous research efforts focus on developing approaches for LWC estimation based on spectral regions from red edge to SWIR [12,13]. There are two widespread methods: (1) the use of spectral indices calculated by the reflectance of wavebands relevant to leaf water status (e.g., water index and normalized difference water index) [14,15], and (2) radiative transfer modeling, that simulates the light penetration path through the canopy, based on physical laws [16]. However, in both these methods it is generally assumed that plant canopies are vertically homogeneous. A major problem arises from the fact that a non-uniformed vertical distribution of many leaf biochemical parameters, including LWC, among different leaf positions has been reported for various crops [17][18][19][20]. This would be one of the principal factors controlling the total light reflected from the canopy [21] and influencing the sensitivity of spectral indices to leaf biochemical parameters. Therefore, particular emphasis needs to be placed on the vertical distribution of LWC when studying on the quantitative models for remote precision monitoring.
Few studies have been carried out to detect the vertical distribution of leaf parameters within canopy using remotely sensed data. The majority of spectral measurements are collected using ground-based spectroscopy, since the sensor is capable of viewing an area of interest with high spectral and spatial resolutions, providing accurate information on the underlying mechanisms of vertical leaf compositions monitoring. In this context, researchers tried to develop estimation models for the vertical distribution of leaf variables [20,22]. Li et al. [23] investigated the contributions of different vertical layers to canopy reflectance spectra by removing wheat tissues gradually from top to the bottom, and concluded that the information about leaves in the bottom layer were usually failed to be captured. However, how deep the sensors can sense is partly dependent on the canopy deep of target of interest. Ciganda et al. [24], working on maize, reported that the red-edge chlorophyll index (CI red edge ) could sense the chlorophyll content of the upper seven to nine leaf layers, when employing a hierarchical regression. From the results of Luo et al. [25], leaf nitrogen (N) content of reed in the top three layers can be accurately quantified based on canopy reflectance. The interrelationship of vertical leaf N distribution within the canopy was simulated by a statistical method and used for the estimation of leaf N content for the whole canopy. In addition, a physically based multiple-layer canopy reflectance model was proposed by Wang et al. [21], and was successfully tested in depicting vertical profiles of leaf variables for winter wheat [26]. They suggested that the penetration characteristics and sensitivity of spectral bands used in the spectral indices should also be considered. Based on these results, several wavebands in the NIR and SWIR regions were identified as effective wavelengths for building spectral indices or estimation models for assessing the vertical leaf N distribution [20,22,27].
For the case of LWC, a bell-shaped vertical trend of LWC distribution has been reported in wheat canopies [23,28]. However, there have been few studies concentrating on the challenging issue of detecting vertical distribution of LWC within canopies based on reflectance data and spectral indices derived [28], especially for the canopy after head emergence stage. Studies have demonstrated that the reflectance of crop canopies could be changed by the emergence of crop spikes at the late growth stages [29,30], which would increase the uncertainty for estimation of LWC vertical distribution. Nevertheless, how wheat spikes affect canopy reflectance across the full spectral regions and the assessment of vertical profile of LWC in crop canopies are still poorly understood. Deep knowledge of these mechanisms issues would facilitate to develop a suitable retrieval strategy, and consequently improve the accuracy of estimation of vertical LWC distribution.
The objectives of this study were (1) to investigate the vertical profile of LWC within winter wheat canopies and identify the interrelationship of LWC among the different vertical layers; (2) to analyze the effects of wheat spikes on canopy spectral reflectance and the performance of published water-related spectral indices in quantifying the vertical LWC in wheat; (3) to estimate vertical distribution of LWC within canopies using a newly proposed indirect induction method, then compare its prediction ability with the conventional estimation method.

Field Experiments
The field experiments were conducted at the National Experiment Site for Precision Agriculture (40 • 10.6 N, 116 • 26.3 E), Beijing, China ( Figure 1). This site has a sub-humid continental monsoon climate, with a mean annual rainfall of 507.7 mm and a mean annual temperature of 13.8 • C. Experiment 1 (Exp. 1) was conducted in 2007. Nine winter wheat cultivars were planted with a density of 3 × 10 6 plants ha −1 and a row spacing of 25 cm, which included Lumai21, Jing411, Jingken49, Jing9843, Jingdong12, 9158, 6211, I-93, Laizhou3279. Nitrogen fertilizer as urea was applied at the pre-planting and the stem elongation stages; 345 kg ha −1 compound fertilizer consists of 15% nitrogen, 15% phosphorus and 15% potassium were used prior to sowing. Each cultivar was planted in a plot with an area of 45 × 10.8 m 2 . Experiment 2 (Exp. 2) was conducted in 2017. Two cultivars of Lunxuan167 and Jingdong18 were investigated. Four N fertilization rates were applied for all cultivars: 0 (N0), 150 kg ha −1 (N150), 300 kg ha −1 (N300), and 450 kg ha −1 (N450). Each cultivar was grown in plots of 15 m × 9 m size, on a silty clay loam soil. Canopy spectral measurement and sampling dates were selected during the late growth periods of wheat, i.e., 9 May (Head emergence, Z54), 19 May (Heading, Z59) and 29 May (Milk-filling, Z73) in Exp.1, 9 May (Head emergence, Z54) and 27 May (Milk-filling, Z73) in Exp. 2. It should be noted that, in Exp. 2, in order to focus on investigating the effects of wheat spikes on canopy spectral reflectance and vertical distribution of LWC estimate, for each plot, wheat spikes were removed from the plant after accomplishing spectral measurement for the entire canopy under the premise of ensuring that all other conditions remained unchanged.

Canopy Spectral Reflectance Measurement
An ASD FieldSpec 3 spectrometer (Analytical Spectral Devices, Boulder, CO, USA) was used to measure the canopy spectral reflectance, with a 25 • field-of-view fiber optics. It records the spectral radiance between 350 and 1050 nm with a sampling interval of 1.40 nm and a resolution of 3 nm, and between 1000 and 2500 nm with a sampling interval of 2 nm and a resolution of 10 nm. The measurements were performed under clear sky conditions between 10:00 am and 14:00 pm (Beijing local time).
Within each plot, three subplots with an area of 1 m × 0.6 m (4 rows × 0.6 m long, 0.6 m 2 ) were randomly selected, canopy spectral reflectance was measured at a height of about 1 m above wheat canopies from the nadir direction. Notably, in Exp. 2, after collecting spectral reflectance of entire canopy, all of wheat spikes in the subplot were carefully cut off, and the spectral reflectance of the remaining canopy without spikes was measured. Each spectral measurement of the two experiments was preceded by a dark current measurement and a white reference measurement was taken before and after canopy spectral measurement, using a 99% white Spectralon ® (Labsphere, Inc., North Sutton, NH, USA) reference panel. Ten scans were determined and averaged to obtain the spectral reflectance of the entire canopy or canopy without spikes for each subplot. Three averaged subplots' spectra were used to represent the reflectance of entire canopy or canopy without spikes for each plot.

Vertical Leaf Water Content Distribution Measurement
All samples of wheat in each subplot within the footprint of canopy reflectance acquisitions, were harvested by cutting off the plants at the soil level, and immediately placed in black plastic bags and kept cool until they were brought back to determine LWC in the laboratory. To quantify vertical profiles of LWC within wheat canopy, we divided the plant into three vertical layers according to various functions of leaves, as shown in Figure 2. The top first leaf (i.e., flag leaf) was assigned to the upper-layer, since its important role in contributing to the major part of photosynthesis of canopy and significantly to grain yield and quality at the late growth stages [31], and the top second leaf was assigned to the middle-layer. In addition, the top third and the leaves below were assigned to the bottom-layer, in consideration of the fact that they were proven to be more sensitive to nutrient and water deficiency [32] and leaves below the top third began to senesce during the milk-filling period, mainly consisting of the yellow and withered leaves, which was particularly pronounced for the N0 treatment. For each subplot, leaves belonging to the same vertical layer were cut from the stems and weighed immediately to obtain the fresh weight (FW) using analytical balance. Then they were dried at 105 • C for 30 min in an oven and subsequently dried at 80 • C until constant weight, dry weight (DW) of leaves were recorded after weighing. Leaf water content (LWC) (%) was calculated using the method of weight ratio [33]. The LWC in a given vertical layer of three subplots were finally averaged to represent the LWC in the corresponding layer of the plot.

Construction of New RRD Type of Spectral Indices
Differences in vegetation canopy geometry and illumination angles during the spectral measurements could cause the confounding effects of scattering that are unrelated to leaf biochemical parameters [43]. The theory of multiple signal correction (MSC) was proposed for separating the biochemical light absorption (e.g., LWC) from the physical light scatter for NIR data [44]. For an individual spectrum or sample, the MSC model was expressed as follows: where R ik represents spectral reflectance for sample i at wavelength k; a i and b i represent the additive and multiplicative effects for sample i respectively; R k represents the "standard" value at wavelength k, which is usually expressed by the average value of all samples at wavelength k; e ik represents residual with respect to all other effects in the spectrum that cannot be modeled by a i and b i for sample i at wavelength k. After the transposition, the scatter corrected spectrum R ik, corrected was calculated by: Equations (2) and (3) suggest that the derivation of the relationship between spectral reflectance and leaf biochemical content requires to eliminate the variations due to the additive offset and the multiplicative effect. However, the extract values of scatter coefficients a i and b i are not unique, they depend on different sample datasets of R k used [43].
In this study, the ratio of reflectance difference (RRD) type of indices was explored, which is formulated as: From the formation of RRD, it is obvious that the scatter coefficients a i and b i will be eliminated, thus be insensitive to the sample datasets. Therefore, it can be expressed in Equation (5) as well.
Because of the excellent performance of the published WI and NDWSI indices in estimating the vertical LWC profile over other spectral indices tested (see Section 3.3), inspired, we optimized them to construct four new RRD type of spectral indices from three or four narrow spectral bands, referred to as WI-3 and WI-4, NDWSI-3 and NDWSI-4, which were defined by the following equations: where 850 nm, 900 nm and 970 nm are the existing wavebands in the WI or NDWSI indices, λ3 and λ4 are the third or/and forth wavebands in the four new RRD type of spectral indices. To find the optimum λ3 and λ4 in each index for estimating vertical distribution of LWC, we calculated the spectral reflectance of every band, or all two bands combinations, from 400-2500 nm, to combine with the reflectance of existing bands, respectively, and then correlated with the LWC in the middle-layer. The λ3 in the WI-3 and NDWSI-3, or the λ3 and λ4 in the WI-4 and NDWSI-4 were consequently selected by the values of coefficient of determination (R 2 ). All the calculations were implemented using MATLAB 8.3 (The MathWorks, Inc., Nat-ick, MA, USA). From our result, λ3 = 1350 nm for WI-3, λ3 = 1200 nm for NDWSI-3, λ3 = 825 nm and λ4 = 1013 nm for WI-4 and NDWSI-4, so the four new narrow-band spectral indices were:

Data Analysis
The datasets collected in Exp. 2 were used to explore the influences of wheat spikes on spectral reflectance and the accuracy of vertical LWC estimation, while datasets collected in Exp. 1 and 2 were used for modeling. The values of coefficients of determination (R 2 ) between LWC in vertical layers and published spectral indices, derived from spectral reflectance of the entire canopy with or without spikes, were calculated, which were referred to as R 2 entirecanopy and R 2 canopywithout spikes , respectively. A relative variation rate (R v , %) was used to compare the performance of each spectral index with respect to vertical LWC before and after spike removal, which was formulated as Equation (14). A higher Rv of R 2 suggests a higher degree of variations in spectral indices, i.e., effects of spikes, on vertical LWC estimation.
where R v of R 2 (%) indicates the change of the accuracy of estimation models for LWC in a given vertical layer before and after spike removal. The R v of R 2 (%) > 0 represents an increase in model accuracy after removing spikes, whereas the R v of R 2 (%) < 0 represents a decrease.
Linear regression was used to model the relationship between spectral indices and LWC in different vertical layers using the SPSS 18.0 software (SPSS Inc., Chicago, IL, USA). In this study, we proposed an indirect induction method to establish models for LWC in the upper-and bottom-layer, by considering the effects of wheat spikes and the interrelationship of vertical LWC within canopies. Their estimation model was as follows: where y indicates the LWC in the upper-layer or the bottom-layer, coefficients m and n indicate the slope and intercept of relationship between middle-LWC and upper-LWC, or middle-LWC and bottom-LWC; y Middle-LWC indicates the LWC in the middle-layer, SI indicates the spectral indices used, coefficients a and b indicate the slope and intercept of estimation model for middle-LWC. The validity of models was assessed based on leave-one-out cross-validation approach. Apart from R 2 and root mean square error (RMSE) were employed to test the performance of the models in vertical LWC prediction, the relative error (RE) and Nash-Sutcliffe efficiency (NSE, calculated as Equation (16)) were used to evaluate the predictive ability of models, which were classified into four and three categories respectively by researchers [45,46], as shown in Table 2.
where y i andŷ i are the measured and estimated LWC in a given vertical layers, y i is the mean of the measured LWC in a corresponding layer.

Vertical Variation of LWC within Wheat Canopies
The vertical profiles of LWC within wheat canopies across different growth stages and different N rates are presented in Figure 3a,b. We can observe that the LWC in the three vertical layers were all gradually decreased after the head emergence stage of wheat. However, for each vertical layer, the treatments with higher N rates (i.e., N300 and N450) had higher LWC in comparison to those with lower N rates (i.e., N0 and N150), although the differences between N300 and N450 were not distinctive, indicating that appropriate application of N fertilizer could improve the LWC to a certain extent. However, for all treatments, there was a similar vertical distribution pattern of LWC, showing great vertical heterogeneity within canopies. Leaves in the middle layer had the largest LWC among the three layers of wheat, while the LWC in the upper layer were slightly lower or higher than that in the bottom layer. Moreover, the middle-LWC was strongly and linearly related to the upper-LWC and bottom-LWC, with correlation coefficient of 0.95 and 0.92, respectively (Figure 3c,d).

Effects of Wheat Spikes on Canopy Spectral Reflectance
To inspect the effects of wheat spikes on spectral reflectance, we compared spectral reflectance of the entire canopy and the canopy without spikes over spectral range from 400 to 2500 nm at different growth stages and N rates (Figure 4). Note that the noises caused by the effect of water vapor in the air above the canopy in the spectral reflectance values, ranged from 1360-1370 nm, 1830-1920 nm and 2440-2500 nm, were removed. As shown in Figure 4a,b, wheat spikes had a remarkable effect on canopy spectra. Judging by the response of canopy reflectance to spikes removal in different spectral regions, the largest changes occurred in the NIR and SWIR regions (730-2500 nm), where the reflectance values of canopy without spikes were generally higher than those of the corresponding entire canopy at the same growth stage and N treatment, and little variation was observed in the visible region (400-730 nm). Through carefully examining the reflectance spectra in Figure 4, we were able to detect that spectral reflectance in the 800-1250 nm, 1400-1780 nm and 2000-2350 nm ranges strongly responded to the removal of spikes, which included the sensitive spectral bands utilized to calculate leaf water-related spectral indices in Table 1. In addition, the comparable extents of differences of canopy reflectance before and after spikes removal were shown between the two late growth stages under the same N treatment (Figure 4a), whereas bigger extents of differences appeared among N fertilization treatments (Figure 4b). Compared to the higher N treatments, reflectance of canopy without spikes under the N0 treatment was considerably different from the original reflectance in visible and SWIR (1400-2500 nm) regions. However, among different N rates, the amplitudes of the differences in reflectance before and after spikes removal within NIR high plateau from 730 to 1300 nm were not distinct, except for the N450 with relatively small difference being observed.

Effects of Wheat Spikes on Relationships between Published Spectral Indices and LWC in Vertical Layers
Given that the wheat spikes contributed differently to canopy spectral reflectance over 400 to 2500 nm (Figure 4), they might have different effects on the estimation accuracy of LWC in different vertical layers. The R 2 entire canopy and R 2 canopy without spikes of linear regression analyses for the relationships between published spectral indices and LWC in the upper-, middle-and bottom-layer are summarized in Table 3. In general, for all vertical layers, the values of R 2 entire canopy were all lower than the corresponding R 2 canopy without spikes , except for the CWI and NIDI, suggesting that the presence of wheat spikes negatively affects the accuracy of spectral indices in assessing the vertical LWC variation. This may present a convincing evidence for the decreased accuracy of leaf biochemical estimation models when the data obtained after the heading stage were pooled [47]. To further quantify the attenuation influence of spikes on LWC estimation in each vertical layer, we calculated the Rv between R 2 entire canopy and R 2 canopy without spikes of each spectral index for the three vertical layers. The results are shown in Figure 5. An important information revealed by Figure 5 is that the spikes had significant, but different, effects on the sensitivity of spectral indices to changes in the vertical LWC within canopies, depending on the vertical layers. Specifically, for almost all spectral indices, the Rv of R 2 in the middle layer yielded the smallest values compared with the upper-and bottom-layer, apart from the WBI/NDVI, LWI and NDII, whose Rv of R 2 in the middle layer ranked the second place. These results indicated that the assessment of middle-LWC was less susceptible to the effect of spikes, whereas the estimation of upper-LWC or bottom-LWC was relative strongly influenced.  In addition, we further compared the performances of spectral indices in LWC determination in different vertical layers before spikes removal, based on datasets from Exp.1 and 2. As presented in Table 4, almost all published spectral indices achieved the best R 2 in the relationship with respect to the middle-LWC, whereas the weakest relationships were found in the bottom layer. In consideration of the significant relationships between middle-LWC vs. upper-LWC and middle-LWC vs. bottom-LWC (Figure 3c,d), it was expected that the LWC in the upper and bottom layers could be indirectly assessed by the remote estimation model of middle-LWC at the late stage of wheat, as a result, reducing the variability due to the wheat spikes and improving vertical LWC estimates within canopies.

Estimation of Vertical LWC Distribution of Wheat Using a Method of Indirect Induction
As illustrated in Table 4 and Figure 5, the best behaviors in terms of both sensitivity to LWC changes of the middle layer and less sensitivity to wheat spikes effects, were shown by the NDWSI, WI, WBI and FWBI1 (which were indicated in bold in Table 4). These indices yielded R 2 values larger than 0.56, while Rv of R 2 values less than 30%. Since the formulas and band combinations of the WI, WBI and FWBI1 showed remarkable consistency (Table 1), we used the WI as a representation. Given the fact that both the additive and multiplicative effects exist in canopy reflectance measurements across different wheat cultivars over growth stages, the WI and NDWSI were optimized by adding the third or/and forth wavebands to construct the four new RRD type of indices, expecting to reduce those effects on LWC estimation based on the MSC theory. Reflectance of each waveband or a large number of two-band combinations over the range of 400-2500 nm, were tested as the R λ3 or R λ3 and R λ4 combination in WI-3, NDWSI-3, WI-4 and NDWSI-4 equations (i.e., Equations (6)-(9)), then related to the middle-LWC. The R 2 value was used to evaluate which waveband or band combination should be selected to develop the best performing RRD type of indices ( Figure 6). As shown in Figure 6a,b, in the selection of the third waveband used in the WI-3 and NDWSI-3, the high significance area all appeared in NIR region. It was worth noting that there are two peaks, λ3 = 1350 nm in the WI-3 and λ3 = 1200 nm in the NDWSI-3, achieving the highest R 2 values for determination of middle-LWC.  Figure 7 shows the comparisons of performances of the optimal RRD type of indices and the corresponding published spectral indices in the middle-LWC estimation. The WI-3 ((R 900 − R 1350 )/(R 970 − R 1350 )) and WI-4 ((R 900 − R 825 ) / (R 970 − R 1013 )) outperformed the WI with increase of R 2 values by 10% and 41%, while the NDWSI-3 ((R 850 − R 970 ) / (R 850 + R 970 − 2R 1200 )) and NDWSI-4 (((R 850 − R 825 ) − (R 970 − R 1013 )) / ((R 850 − R 825 ) + (R 970 − R 1013 ))) were superior over the NDWSI with increase of R 2 values by 3% and 38%, respectively. Being more effective than the three-band RRD indices, the two four-band RRD indices reduced the saturation of spectral indices to some degree and increased the accuracy of LWC estimation, where the WI-4 and NDWSI-4 explained 82% and 84% of the variation in middle-LWC at the late stages. Since the significant and linear relationships between middle-LWC vs. upper-LWC and middle-LWC vs. bottom-LWC (Figure 3c,d), along with the distinct sound behaviors of the WI-4 and NDWSI-4 indices which were clearly identifiable in Figure 7 with respect to the middle-LWC, we developed estimation models through a method of indirect induction for the upper-LWC and bottom-LWC respectively as follows: y Bottom−LWC = 0.7447 × y Middle−LWC + 12.276 where y U pper−LWC , y Bottom−LWC and y Middle−LWC indicate the LWC in the upper-, bottomand middle-layer, respectively; SI indicates the WI-4 or NDWSI-4 index, which were formulated as Equations (11) and (13); coefficients a and b indicate the slope and intercept of estimation model based on the WI-4 or NDWSI-4, as shown in Figure 7c,f.

Validation of Vertical LWC Distribution Models
To validate the predictive ability and acceptability of the estimation models of vertical LWC described in Section 3.4, the estimated LWC in the middle-, upper-and bottom-layer, derived from models based on the WI-4, NDWSI-4 and the best performing published ND-WSI, were compared to the measured LWC in the corresponding vertical layers (Figure 8). It can be observed that for all plots, all REs were less than 30%, NSEs were larger than 0.5, and almost all scatters fell into the 95% confidence intervals of prediction, confirming the acceptability of all models tested. The WI-4 and NDWSI-4 predictions for the LWC in the three vertical layers exhibited advantages over the published NDWSI (R 2 < 0.6), which had led to very higher coefficients of determination (R 2 ≥ 0.74). Specifically, for the middle-layer, we have found a very consistent agreement between LWC values measured and those estimated by WI-4 and NDWSI-4, with R 2 and NSE values higher than 0.8, and REs lower than 10%. The NDWSI-4 models performed the best in estimating the LWC in the top two layers, which explained 83% and 81% the variations in middle-LWC and upper-LWC, respectively. For the bottom layer, the WI-4 and NDWSI-4 produced similar result, with R 2 of 0.75 and 0.74. However, the NDWSI-4 estimate was more preferable, since it generated lower RMSE and RE values, as well as larger NSE, expressing distribution of scattering points being closer to the one-to-one line than the WI-4 estimate. In order to evaluate whether the indirect induction method proposed in this study has merit in vertical LWC prediction, we compared it with the conventional direct estimation method. Linear estimation models between the upper-LWC and WI-4 and NDWSI-4, as well as between the bottom-LWC and WI-4 and NDWSI-4, were directly developed without taking into account the spikes effects and internal relationship of LWC among vertical layers. The model validations based on the two methods for the upper and bottom layers are shown in Figure 9. The most encouraging result is that, for the two layers, all LWC models established by the indirect induction method outperformed the models developed by the direct estimation, indicated by the increase in R 2 from 0.66 to 0.78 and from 0.77 to 0.81 for the models based on the WI-4 and NDWSI-4 for the upper-LWC, while the increase in R 2 from 0.68 to 0.75 and from 0.67 to 0.74 for the bottom-LWC, respectively.

Discussion
After the head emergence stage of winter wheat, wheat spikes grow out and the canopies become more and more dense, literature shows that leaf biochemical variables based on remote sensing technique perform well in the early stage, however, it tends to produce uncertainty or error during the reproductive stage [48], and then could result in more difficulties to monitor their vertical profile within canopies. Thus, it is vital to fully understand the effects of spikes on canopy reflectance and the detection degree of leaf information in different vertical layers by remotely sensed data, then to derive spectral indices to estimate vertical leaf biochemical parameters. In this work, we explored the remote estimation of the vertical distribution of LWC, considering the effects of spikes, using field data of winter wheat. From the chemical analyses determined in the laboratory, the variation trend in LWC, with higher values in the middle leaves than in the upper and bottom leaves across different growth phases and N rates, was in agreement with the results of Li et al. [23]. This is mostly because water is transported from the bottom shaded leaves to the top sunlit leaves, in an effort to provide materials to reach the maximum total canopy photosynthesis, but the top leaves (i.e., the fully expanded flag leaves) exhibited less LWC due to being more directly exposed to solar irradiance and consume more water through evapotranspiration [28].
The effects of spikes and other plant components (e.g., stem) on canopy spectral characteristics have gradually attracted attention, relevant researches have been conducted using wavelengths lower than 1000 nm [30,49]. However, spectral discrepancies across the whole spectrum required concern for an improved understanding of the interaction between spectral reflectance and leaf biochemical parameters across the canopy vertical profile. In this study, we investigated canopy spectral characteristics, before and after spikes removal, from 400 nm to 2500 nm, across different N rates, wheat cultivars and growth stages. Our results showed that the values of reflectance for the entire canopy were largely lower than those for the canopy without spikes at NIR and SWIR regions, which was similar to the findings of Haboudane et al. [8] and He et al. [50]. Based on the theory reported in literature, that a beam of radiation becomes diffuse as it passes through the crop canopy, and when it encounters different geometric structures of plant tissues, it will be scattered at each refractive discontinuity [51], the results from our study could be mainly explained by the fact that the emergence of spikes changed the original canopy structures, resulting in changes of the radiation path and light distribution within wheat canopies, and consequently influenced the multiple scattering effects of canopy reflectance, thereby inducing a decrease in spectral reflectance.
In comparison to the visible region, the emerged spikes had larger effects on the canopy reflectance ranged from NIR to SWIR (730-2500 nm), expressing larger variations in these regions (Figure 4). It should be noted that the reflectance of wavebands sensitive to leaf water or water stress were included as well. These bands were usually employed in leaf water-related spectral indices, in turn, led to the significant differences in their performances in estimating the vertical LWC before and after removing spikes. Gutierrez et al. [30] and He et al. [50] found that the spike effects could alter and make weaker relationships between spectral indices and crop variables (e.g., crop yield and LAI), our results showed similarities and differences with such findings through the controlled experiment for either remaining or removing wheat spikes. The appearance of spikes decreased the accuracy of the vertical distribution of LWC estimation, but there were differences in the degrees of the effects depending on the vertical layers, since the superior insensitivity of middle-LWC estimation to the spike effects were revealed, compared to the upper-and bottom-LWC estimation (especially for the NDWSI, WI, WBI and FWBI1 in Table 3). In addition, better correlation between spectral indices and the middle-LWC than the other two layers shown in this work (Table 4) were consistent with the results of Liu et al. [28], this provides an evident for differences of ability of the nadir remote sensing data in detecting vertical LWC. By comparing the finding reported in Li et al. [32], that the reflectance of leaves positioned in the upper layer contributed relative more to the reflected light gathered by sensors, the result of this study is probably because the spikes negatively affected the LWC information acquisition of upper leaves after head emergence stage, thus tending to obtain lower estimation accuracy than the middle-LWC.
Compared with the estimation of crop biochemical variables at the leaf or canopy scales, the estimation of their vertical distribution within canopies is challenging, especially during growth periods after spikes appearance. Some studies have attempted to explore a multi-angular remote sensing technology and conduct a great deal of work for enhancing the vertical variables retrieval [19,52]. Although this approach can provide more information on canopy vertical structure by viewing from different angles than that only from a single nadir direction, there are still shortcomings. For instance, complex operation and time-consuming operation of multi-angle spectral measurements [53], as well as the added uncertainty resulting from view angle effects on spectral indices [54]. Jacquemoud et al. [55] reported that the accuracy of vegetation parameters estimation depends largely on the retrieval technique. In this study, we tried to detect the vertical distribution of LWC within wheat canopies based on canopy reflectance obtained from the nadir observation, using a new method of indirect induction. Firstly, by consideration of the efficacy of remote sensing in detecting vertical LWC and the wheat spikes effects on the estimation of LWC in each vertical layer, we identified the most effective leaf layer (i.e., the middle layer) in which the spectral indices not only achieved the highest sensitive to the corresponding LWC but also less susceptible to the emergence of spikes, and developed the estimation model for this layer based on the MSC theory. Subsequently, accompanied with the significant internal relationship of LWC among vertical layers, the LWC in the remaining layers were indirectly deduced by the estimation model of the above effective layer. In the vertical LWC determination, our optimized three-or four-narrow-band spectral indices, WI-3, WI-4, NDWSI-3 and NDWSI-4, presented substantially better behavior in detecting the LWC across all the treatments than the original indices, in particular for the WI-4 and NDWSI-4 (R 2 higher than 0.81 in middle-LWC assessment). These four new spectral indices were designed by adding one or two NIR bands into the WI and NDWSI on the basis of RRD formulation. Their advantage probably relies on the compensation of additive and multiplicative effects resulted from different canopy structures, anisotropic multiply scattering and soil background by using the MSC method, because the effects are common in the measurements of canopy spectral reflectance but have nothing to do with leaf biochemical parameters, i.e., LWC in this study. Similarly, the formula of RRD was also used by researchers to develop spectral indices, such as the modified SR (mSR 705 = (R 750 − R 445 )/(R 705 − R 445 )) and modified NDVI (mND 705 = (R 750 − R 705 )/(R 750 + R 705 − 2R 445 )), the MERIS terrestrial chlorophyll index (MTCI = (R 754 − R 709 )/(R 709 − R 681 )), the structure-insensitive pigment index (SIPI = (R 800 − R 445 )/(R 800 − R 680 )) and the modified difference ratio (MDR = (R 1271 − R 410 )/(R 1342 − R 410 )) [7,[56][57][58], they all have been proven to be an improvement in terms of sensitivity to leaf pigment and leaf water contents when applied across a wide range of plant species, leaf structures and growth stages. From the comparison results between methods of the indirect induction and the conventional direct estimation (Figure 9), we concluded that the newly proposed approach contributed to achieve the higher potential of canopy spectra obtained from the nadir direction in monitoring vertical profiles of LWC after head emergence of wheat. Due to the limitation of samples, this method needs to be further studied by more experimental data and crop cultivars with different geometry types in the future. Although an increasing number of optical satellite images are freely available, the use of space-based imaging spectroscopy to quantify the vertical LWC distribution within crop canopies is still relatively challenging, because of the limited spatial and spectral resolutions. However, the application of narrow-band spectral indices to detect leaf biochemical parameter details has been an intended goal for studies of canopy physiology and ecology, our results provides support for the selection of spectral bands to design future sensors, and eventually to advance the identification and mapping of the vertical LWC from satellite data.

Conclusions
This study shows that leaf water content (LWC) within winter wheat canopies tended to be higher in the middle layer and lower in the upper and bottom layers, the middle-LWC yielded very strong correlations with the upper-LWC and bottom-LWC after head emergence stage. The effects of wheat spikes on canopy reflectance were analyzed. We found that the presence of spikes had significant effects on the reflectance in NIR to SWIR spectral regions, with the water absorption bands included. In addition, the spikes effects negatively impacted on the estimation of LWC in the vertical layers. However, LWC in the middle layer was identified to be the most highly sensitive to almost all earlier proposed water-related spectral indices and be more resistant to wheat spikes effects, showing an advantage over the upper and bottom layers. This implied the potential of applying middle-LWC estimation models to extend the determination of LWC also in the other layers. Based on the results, a method of indirect induction was proposed to establish models of vertical distribution of LWC, by considering the effects of wheat spikes. The narrow-band WI-4 ((R 900 − R 825 ) / (R 970 − R 1013 )) and NDWSI-4 (((R 850 − R 825 ) − (R 970 − R 1013 )) / ((R 850 − R 825 ) + (R 970 − R 1013 ))) indices were developed based on the multiply signal correction theory and achieved the highest sensitivity in capturing variations in the middle-LWC. Importantly, the prediction accuracy of LWC in the upper and bottom layers was shown to be improved by the indirect models derived from these two indices, compared to the direct estimation. Future studies on more crop or vegetation species will be of importance for validating and improving the method developed here.
Author Contributions: W.K. processed and analyzed the field measurements and hyperspectral data, and wrote the manuscript. W.H., L.T. and C.L. guided the data analysis and provided suggestions for the study. L.M. and X.Z. designed the experiment and were involved in the ground data collection. R.C. revised, provided suggestions and edited the manuscript. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The datasets referred to in this paper are available from the corresponding author on reasonable request.