Assessment of Leaf Area Index Models Using Harmonized Landsat and Sentinel-2 Surface Reflectance Data over a Semi-Arid Irrigated Landscape

Leaf area index (LAI) is an essential indicator of crop development and growth. For many agricultural applications, satellite-based LAI estimates at the farm-level often require near-daily imagery at medium to high spatial resolution. The combination of data from different ongoing satellite missions, Sentinel 2 (ESA) and Landsat 8 (NASA), provides this opportunity. In this study, we evaluated the leaf area index generated from three methods, namely, existing vegetation index (VI) relationships applied to Harmonized Landsat-8 and Sentinel-2 (HLS) surface reflectance produced by NASA, the SNAP biophysical model, and the THEIA L2A surface reflectance products from Sentinel-2. The intercomparison was conducted over the agricultural scheme in Bekaa (Lebanon) using a large set of in-field LAIs and other biophysical measurements collected in a wide variety of canopy structures during the 2018 and 2019 growing seasons. The major studied crops include herbs (e.g., cannabis: Cannabis sativa, mint: Mentha, and others), potato (Solanum tuberosum), and vegetables (e.g., bean: Phaseolus vulgaris, cabbage: Brassica oleracea, carrot: Daucus carota subsp. sativus, and others). Additionally, crop-specific height and above-ground biomass relationships with LAIs were investigated. Results show that of the empirical VI relationships tested, the EVI2-based HLS models statistically performed the best, specifically, the LAI models originally developed for wheat (RMSE:1.27), maize (RMSE:1.34), and row crops (RMSE:1.38). LAI derived through European Space Agency’s (ESA) Sentinel Application Platform (SNAP) biophysical processor underestimated LAI and provided less accurate estimates (RMSE of 1.72). Additionally, the S2 SeLI LAI algorithm (from SNAP biophysical processor) produced an acceptable accuracy level compared to HLS-EVI2 models (RMSE of 1.38) but with significant underestimation at high LAI values. Our findings show that the LAI-VI relationship, in general, is crop-specific with both linear and non-linear regression forms. Among the examined indices, EVI2 outperformed other vegetation indices when all crops were combined, and therefore it can be identified as an index that is best suited for a unified algorithm for crops in semi-arid irrigated regions with heterogeneous landscapes. Furthermore, our analysis shows that the observed height-LAI relationship is crop-specific and essentially linear with an R2 value of 0.82 for potato, 0.79 for wheat, and 0.50 for both cannabis and tobacco. The ability of the linear regression to estimate the fresh and dry above-ground biomass of potato from both observed height and LAI was reasonable, yielding R2: ~0.60.


Study Area and Surveyed Crops
The study area is focused on the agricultural scheme of the Bekaa Valley of Lebanon, which is an almost flat plain between two mountain ranges ( Figure 1). The total cropland area in the valley is 118,000 ha [45], which accounts for 42% of the total cultivated area in Lebanon [46,47]. The agricultural sector in Lebanon is a key contributor to the country's agri-food industry. Most of the rural population in the underprivileged regions of Lebanon (Akkar, Dinniyeh, the Northern Bekaa, and the South) rely on agriculture activities, which account for 80% of the local GDP [48]. The study area covers Baalbek (characterized by semi-arid climate) in the 2019 growing season and the North to West Bekaa (characterized by the Mediterranean climate) in the 2018 growing season. Field campaigns were conducted within the study area to acquire ground LAI measurements in several different crop types. In total, 451 fields were surveyed representing a total cultivated area of 5606 ha. A summary of the fields and crops surveyed and measurements collected over the two growing seasons is given in Table 1. Of the crops studied here, wheat and potato are the two primary strategic cash crops of the Bekaa Valley. In addition, vegetables cover around 20% of the agricultural lands in Lebanon [49]. Due to the involvement of large numbers of people in the cultivation, processing, and trading of tobacco, this crop plays a vital role in the economic and financial sector of the country. Of the total tobacco produced, 43% is produced in the Bekaa and the north of Lebanon [50]. Additionally, cannabis is widely cultivated in Central and Northern Bekaa, with more than 7000 ha planted annually [51]. Cannabis has a high return value compared to other crops representing a major source of income for some farmers in the Bekaa Valley.

Leaf Area Index Measurements
In situ measurements of LAI were collected in the 2018 and 2019 growing seasons using the SS1 SunScan Canopy Analysis System (Delta-T Devices). This instrument consists of 64 sensors that measure photosynthetic active radiation (PAR) with an accuracy of ±10%, embedded in a 1 m long portable probe (Figure 2a) with an external beam fraction sensor (BFS) (Figure 2b) [52]. The BFS is used to monitor the light incident on the canopy at the same time as the below-canopy measurements are taken. The SS1 requires several inputs, including the ellipsoidal leaf angle distribution (ELADP), longitude, latitude, and local timing (Figure 2a). Typical ELADP values used in this study are summarized in Table A1 in Appendix A. LAI measurements were taken during the same day as a Landsat 8 or Sentinel 2 satellite overpass to calibrate and evaluate the relationship between LAI and the vegetation indices to be tested. Some measurements were taken within ±2 h of the overpass, and others continued to the rest of the day.  Two different patterns for data collection strategies were used in 2018 and 2019. In 2018, LAI measurements were collected over the full Bekaa Valley study area in Figure 1. To capture withinfield variability, five measurements were taken in each field, and the average value was taken to be representative of the field-scale. For each field, the measurements were acquired at five different  Two different patterns for data collection strategies were used in 2018 and 2019. In 2018, LAI measurements were collected over the full Bekaa Valley study area in Figure 1. To capture within-field variability, five measurements were taken in each field, and the average value was taken to be representative of the field-scale. For each field, the measurements were acquired at five different locations~5 m apart following the pattern schematized in Figure 3 (points A, B, C, D, and E), starting from the center of the field. Measurements were collected away from the field boundaries to avoid edge effects. Two different patterns for data collection strategies were used in 2018 and 2019. In 2018, LAI measurements were collected over the full Bekaa Valley study area in Figure 1. To capture withinfield variability, five measurements were taken in each field, and the average value was taken to be representative of the field-scale. For each field, the measurements were acquired at five different locations ~5 m apart following the pattern schematized in Figure 3 (points A, B, C, D, and E), starting from the center of the field. Measurements were collected away from the field boundaries to avoid edge effects. For the 2019 growing season, extending from April to July, the sampling focus was on monitoring time-series LAI dynamics of three significant crops in the Bekaa regions: wheat, barley, and potato. A total of 42 predefined sampling sites in North Bekaa with recorded GPS coordinates were selected, and LAI was measured at these locations on each Landsat 8 or S2 overpass date. The sampling sites include 21 locations in wheat fields, 12 locations in barley fields, and 9 locations in potato fields.

Canopy Height and Above-Ground Dry and Fresh Biomass Measurements
Measurements of above-ground crop biomass and height were also collected at a subset of the LAI sampling sites in 2018 and 2019. A quadrate (0.5 m 2 × 0.5 m 2 ) was randomly thrown in the field, For the 2019 growing season, extending from April to July, the sampling focus was on monitoring time-series LAI dynamics of three significant crops in the Bekaa regions: wheat, barley, and potato. A total of 42 predefined sampling sites in North Bekaa with recorded GPS coordinates were selected, and LAI was measured at these locations on each Landsat 8 or S2 overpass date. The sampling sites include 21 locations in wheat fields, 12 locations in barley fields, and 9 locations in potato fields.

Canopy Height and Above-Ground Dry and Fresh Biomass Measurements
Measurements of above-ground crop biomass and height were also collected at a subset of the LAI sampling sites in 2018 and 2019. A quadrate (0.5 m 2 × 0.5 m 2 ) was randomly thrown in the field, and the vegetation within the frame of the quadrate was clipped. Three samples of potato biomass were taken from each field, and the obtained results were averaged. The above-ground plant component was oven-dried at 75 • C for 48 h to a constant weight and then weighed. Above-ground biomass measurements were collected in the potato fields in 2019. Crop height was directly measured for the tobacco and cannabis crops in 2018, and the potato crop in 2019 for the plant samples within the frame. Figure 1 shows the biophysical parameters sampling sites (labeled in dark green) for the 2019 growing season. Crop height and LAI field sampling at different growing stages of a wheat culture studied during the 2019 growing season is shown in Figure 4.  Figure 1 shows the biophysical parameters sampling sites (labeled in dark green) for the 2019 growing season. Crop height and LAI field sampling at different growing stages of a wheat culture studied during the 2019 growing season is shown in Figure 4.

Satellite Datasets
Satellite data were obtained for dates when measurements were collected in the field during the 2018 and 2019 growing seasons. Additional scenes were collected in 2019 to develop the time-series of biophysical retrievals. A total of 17 and 24 HLS images (Tile T36 SYC) were used for the 2018 and 2019 growing seasons, respectively. HLS has adopted the tiling system used by Sentinel-2. The tiles are in the Universal Transverse Mercator (UTM) projection and are nominally 110 km on a side. The tiling system is aligned with the UTM-based Military Grid Reference System (MGRS). The first two numbers of a tile name (such as 36SYC) correspond to the UTM zone, 6° of longitude in width. Each zone is divided into latitude bands of 8°. This is represented by a letter, S in the tile covering the study area. Finally, each 6° × 8° polygon (grid zone) is further divided into 110 km tiles, from west to east,

Satellite Datasets
Satellite data were obtained for dates when measurements were collected in the field during the 2018 and 2019 growing seasons. Additional scenes were collected in 2019 to develop the time-series of biophysical retrievals. A total of 17 and 24 HLS images (Tile T36 SYC) were used for the 2018 and 2019 growing seasons, respectively. HLS has adopted the tiling system used by Sentinel-2. The tiles are in the Universal Transverse Mercator (UTM) projection and are nominally 110 km on a side. The tiling system is aligned with the UTM-based Military Grid Reference System (MGRS). The first two numbers of a tile name (such as 36SYC) correspond to the UTM zone, 6 • of longitude in width. Each zone is divided into latitude bands of 8 • . This is represented by a letter, S in the tile covering the study area. Finally, each 6 • × 8 • polygon (grid zone) is further divided into 110 km tiles, from west to east, 4th letter (Y in this case), and south to north, 5th letter (C in this case). For further information on the output projection and gridding, readers are referred to the HLS product user's guide [39]. Figure 5 represents a summary of the satellite data used in the study.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 39 output projection and gridding, readers are referred to the HLS product user's guide [39]. Figure 5 represents a summary of the satellite data used in the study.  Imagery from the newly released harmonized surface reflectance product (HLS, v1.4) was downloaded for these dates from May to September in 2018 and March to August in 2019. HLS represents a synergy of Landsat-8 (L30) and Sentinel-2 (S30) data with 30 meters resolution [39]. This harmonized data product is based on the framework established by Claverie et al. [39], which encompasses several steps: resampling Sentinel-2 to 30 m, a bidirectional reflectance distribution function (BRDF) normalization for both Landsat 8 and Sentinel 2 datasets, a geographical registration for the two products, and a bandpass adjustment using the Landsat-8 band passes as a reference to which the Sentinel-2 spectral bands are adjusted. The resulting data products, S-30 MSI and L-30 OLI consist of surface reflectance data that are consistent in terms of spectral response and spatial resolution. A detailed representation of the HLS processing methodology can be found in [39].
Atmospherically corrected Level 2-A surface reflectance images from Sentinel-2 were also downloaded from the THEIA website at the French Land Data Center (https://www.theia-land.fr/) Imagery from the newly released harmonized surface reflectance product (HLS, v1.4) was downloaded for these dates from May to September in 2018 and March to August in 2019. HLS represents a synergy of Landsat-8 (L30) and Sentinel-2 (S30) data with 30 m resolution [39]. This harmonized data product is based on the framework established by Claverie et al. [39], which encompasses several steps: resampling Sentinel-2 to 30 m, a bidirectional reflectance distribution function (BRDF) normalization for both Landsat 8 and Sentinel 2 datasets, a geographical registration for the two products, and a bandpass Remote Sens. 2020, 12, 3121 8 of 35 adjustment using the Landsat-8 band passes as a reference to which the Sentinel-2 spectral bands are adjusted. The resulting data products, S-30 MSI and L-30 OLI consist of surface reflectance data that are consistent in terms of spectral response and spatial resolution. A detailed representation of the HLS processing methodology can be found in [39].
Atmospherically corrected Level 2-A surface reflectance images from Sentinel-2 were also downloaded from the THEIA website at the French Land Data Center (https://www.theia-land.fr/) on the S30 dates listed in Figure 5. These were used to generate LAI with the Sentinel LAI index and SNAP software, designed specifically for use with THEIA SR data products. The goal was to use the LAI biophysical product from Sentinel-2 (using SNAP and S2 index) as a benchmark to compare to the LAI generated using HLS.

Cloud Masking
Cloud masking is required to identify high-quality land surface reflectance values with low uncertainties. The HLS product contains, in addition to the surface reflectance bands, a level-2-pixel quality assurance band. This band provides information about the cloud status in each pixel (clear, shadow, low, medium, and high confidence), snow/ice, and water pixels. The Landsat-8 cloud mask is derived from both the mask in the USGS Landsat TOA data and the Land Collection 1 Surface Reflectance Code (LaSRC) atmospheric correction tool. However, the Sentinel-2 cloud mask is a union of the Fmask algorithm, which has been adapted from [53] and the LaSRC mask. The quality assessment band (QA) bits of the HLS Landsat and Sentinel-2 were decoded with a loop using simple integer arithmetic, and its output was used to mask cloud shadow, adjacent cloud, cloud, and cirrus from the processed images.
As for the L2A products from THEIA, atmospheric correction, and detection of clouds and cloud shadows are performed using the Multi-sensor Atmospheric Correction and Cloud Screening (MACCS) processor (https://labo.obs-mip.fr/multitemp/?p=6203). The L2A 10-m cloud mask product was downloaded and used from the mask directory.

Vegetation Indices
Vegetation indices (VIs) are widely used in the assessment of several biophysical parameters including the fraction of photosynthetically active radiation absorbed by the vegetation (FAPAR) [54], green vegetation fraction [55], leaf area index (LAI) [56], and canopy chlorophyll content [57]. In this study, we use several indices for the evaluation of LAI, including the soil adjusted vegetation index (SAVI), normalized difference vegetation index (NDVI), enhanced vegetation index 2 (EVI2), and the Sentinel-2 LAI Index (SeLI). SAVI, NDVI, and EVI2 are a combination of visible and near-infrared bands, while SeLI is a combination of near-infrared and the red edge bands of the Sentinel 2 imagery.
NDVI was selected as one of the simplest and earliest VI. Compared to NDVI, both SAVI and EVI2 are less sensitive to soil background [58], and EVI2 is less saturated at high LAI ranges [58]. For SAVI, we use L = 0.5 (USGS, 2017b). SeLI is a new robust LAI green index developed by Pasqualotto et al. [30] and computed from near-infrared and red edge bands. The red-edge region, between 690 and 750 nm, represents the region between the maximum absorption in the red wavelength and maximum reflection in the near-infrared wavelength by the plant [59,60]. Several researchers emphasize the importance of the red edge band in the estimation of LAI green . Equations for the VIs used in this study are presented in Table 2, along with primary citations.
In the literature, vegetation indices are often calculated from top-of-atmosphere reflectance. In this study, however, we calculate the VIs (NDVI, SAVI, and EVI2) using the HLS (Harmonized Landsat Sentinel-2) 30-m surface reflectance products derived from both Landsat-8 L1T products and Sentinel-2 MSI L1C data. Furthermore, we calculated SeLI from the L2A THEIA surface reflectance product. Several existing VI-based LAI models were used to create spatial LAI layers for the study area (Table 3). These models were evaluated with the collected LAI measurements to determine the formula with the best performance over the sampled crop types. The models evaluated include relationships with NDVI (linear, exponential, and polynomial of third-order fits), SAVI (third-order polynomial and logarithmic fits), EVI2 (square root fits), and linear fits to SeLI. The relationships between LAI and different studied vegetation indices, as presented in literature, are summarized in Table 3.

SNAP Model
In addition to the VI-based LAI models described in Section 2.4.1, LAI was produced from S2 THEIA surface reflectance data using Sentinel Application Platform toolbox (SNAP) software. SNAP provides a "Biophysical Processor" tool in the thematic land processing pull-down menu, used for the retrieval of several biophysical variables, including LAI. The algorithm is based on neural networks that are trained to estimate the canopy characteristics [66]. LAI can be directly retrieved for each pixel based on a pre-trained neural net. The neural nets are trained using well-known RTM (PROSAIL: PROSPECT [67] + SAIL [26]) as described in [28]. The ANN algorithm requires eight S2 spectral wavebands (B3-B7, B8a, B11, and B12) as input. The eight bands are resampled to 20 m to preserve the red-edge region in the S-2 scenes. LAI, with the quality indicators of the used dataset, were thus automatically obtained with the finest possible resolution. Figure 6 represents a summary of the overall methodological approach that was followed in this study to generate LAI from the studied vegetation indices computed from HLS and the THEIA L2A surface reflectance products and using the SNAP software.  [26]) as described in [28]. The ANN algorithm requires eight S2 spectral wavebands (B3-B7, B8a, B11, and B12) as input. The eight bands are resampled to 20 m to preserve the red-edge region in the S-2 scenes. LAI, with the quality indicators of the used dataset, were thus automatically obtained with the finest possible resolution. Figure 6 represents a summary of the overall methodological approach that was followed in this study to generate LAI from the studied vegetation indices computed from HLS and the THEIA L2A surface reflectance products and using the SNAP software.

Figure 6.
Step-by-step flow chart for the validation of satellite-derived LAI against the ground measured LAI. Components of the analysis include: (1) existing vegetation index (VI) relationships applied to Harmonized Landsat-8 and Sentinel-2 (HLS) surface reflectance generated by NASA, (2) THEIA L2A surface reflectance products from Sentinel-2, (3) the SNAP biophysical model, (4) groundtruthing, (5) statistical comparison between the measured and the modeled LAI, (6) monitoring crop biophysical parameters during the 2019 growing season, and (7) investigation of crop-specific height and above-ground biomass relationships with LAI.
In this methodological approach, we first perform a comparative assessment of LAI models (both physical and existing empirical LAI estimation methods) against the ground LAI measurements of the studied crop groups. Particular focus is placed on evaluating the newly embedded biophysical processor, based on the artificial neural network (ANN) method, in the Sentinel Application Platform (SNAP) software, as well as the newly developed S2 LAI algorithm from SeLI index, in comparison with empirically-based VI models applied to HLS data. Next, we generated new empirical fits of the VI data to our ground observations specific to these crop groups. Finally, we see how well LAI (observed and retrieved) correlates with other important crop development variables: canopy height and biomass. To further assess the performance of the studied LAI models against the measured LAI, the measured LAI data were divided into three groups: low (0-2), medium (2)(3)(4), and high (>4) LAI ranges.

Figure 6.
Step-by-step flow chart for the validation of satellite-derived LAI against the ground measured LAI. Components of the analysis include: (1) existing vegetation index (VI) relationships applied to Harmonized Landsat-8 and Sentinel-2 (HLS) surface reflectance generated by NASA, (2) THEIA L2A surface reflectance products from Sentinel-2, (3) the SNAP biophysical model, (4) ground-truthing, (5) statistical comparison between the measured and the modeled LAI, (6) monitoring crop biophysical parameters during the 2019 growing season, and (7) investigation of crop-specific height and above-ground biomass relationships with LAI.
In this methodological approach, we first perform a comparative assessment of LAI models (both physical and existing empirical LAI estimation methods) against the ground LAI measurements of the studied crop groups. Particular focus is placed on evaluating the newly embedded biophysical processor, based on the artificial neural network (ANN) method, in the Sentinel Application Platform (SNAP) software, as well as the newly developed S2 LAI algorithm from SeLI index, in comparison with empirically-based VI models applied to HLS data. Next, we generated new empirical fits of the VI data to our ground observations specific to these crop groups. Finally, we see how well LAI (observed and retrieved) correlates with other important crop development variables: canopy height and biomass. To further assess the performance of the studied LAI models against the measured LAI, the measured LAI data were divided into three groups: low (0-2), medium (2)(3)(4), and high (>4) LAI ranges.

Accuracy Assessment
Statistical analysis between observed and remotely sensed LAI was performed for each crop group during the 2018-2019 growing seasons. In 2018, the study was conducted for five crop groups, namely: herbs, grains, potato, tobacco, and vegetables. In 2019, the study was conducted on three major crops, namely: wheat, barley, and potato.

of 35
Root mean square error (RMSE) [68]: Mean absolute error (MAE) [69]: Mean absolute percentage error (MAPE) [70]: Relative RMSE over the average measured LAI, also referred to as coefficient of variation of root-mean square error-CV(RMSE) as in [71]: Mean bias error (MBE) [72]: Index of agreement (d) [73]: We also investigate the relation between the measured LAI and the remotely sensed vegetation indices. For this purpose, we use Pearson's r, an often-used parametric correlation coefficient given by [74]: The significance of the correlation is computed using the two-tailed t-test with n-2 degrees of freedom:

Results
The performance of the empirical and biophysical models was evaluated using the aforementioned statistical parameters in comparison with field measurements. A table of the derived statistical metrics, along with a summary graphical representation, is shown in Figure 7 and Table 4, respectively. The results are explained in the sections below. , and the region between these two vertical lines represents the medium observed LAI range (2-4 m 2 /m 2 ), the 2-D density plots represent four color levels with a dark red color at high density to dark green at lower densities, in addition to the orange and light green in representing medium densities. , and the region between these two vertical lines represents the medium observed LAI range (2-4 m 2 /m 2 ), the 2-D density plots represent four color levels with a dark red color at high density to dark green at lower densities, in addition to the orange and light green in representing medium densities.

HLS VI Models
Empirical LAI models that are developed in relation to vegetation indices are considered as one of the most straightforward methods used for predicting LAI. Here we evaluate the efficiency and robustness of LAI models developed from indices that are computed from red and near-infrared HLS surface reflectance bands using the normalized difference ratio: NDVI (Model 1, 2, and 3), soil adjusted vegetation index: SAVI (Model 4 and 5), and enhanced vegetation index 2: EVI2 (Model 6, 7, 8, and 9). Our analysis (Figure 7 and Table 4) shows that Model 7 (EVI2, initially developed for row crops), Model 8 (EVI2, developed for maize), and Model 9 (EVI2, developed for wheat) outperformed other models when LAI observations from all crops sampled in the current study are combined.
On the other hand, models performing relatively poorly can also be identified. Both Model 2 (NDVI, initially developed for wheat, and corn) and Model 4 (SAVI, for overall agricultural crops) underestimated the LAI for all crops. Underestimation was also noticed in Model 5 (SAVI, overall agricultural crops) and Model 6 (EVI2, overall agricultural crops) for grains. Model 8 (EVI2, initially developed for maize) and Model 3 (NDVI, overall agricultural crops) overestimated the LAI for potato. Our analysis shows an apparent underestimation of LAI in all tested models when measured LAI exceeded 4, except for Model 3 and Model 8 which overestimated LAI for potato for observed LAI > 4.

Sentinel-2 LAI Models
For comparison with SNAP and SeLI LAI derived from S2 data only, subsets of HLS model-observation pairs were selected on the S2 dates represented in the THEIA dataset ( Figure 8 and Table 5). The highest performing VI models applied to HLS SR data (7, 8, and 9; Section 3.1.1) outperformed LAI derived from S2 SeLI and the SNAP model (see Section 3.1.2; Table 5 and Figure 8). The HLS VI models had lower RMSE, MAE, RRMSE, and MAPE values than those of Model 10 (SeLI, overall crops), and the SNAP model. Model 9 (EVI2, initially developed for wheat) had the highest index of agreement (d), lowest RMSE, MAE, and MBE, compared to other models, primarily Model 10 (SeLI, overall crops) and the SNAP model. Furthermore, LAI derived from HLS EVI2 outperformed the LAI models derived from the other studied vegetation indices.

Empirical Relationships between VIs and LAI Observations in Bekaa
Vegetation indices are shown to correlate well with LAI. The studied VIs were evaluated with the multi-crop observed LAI dataset over both 2018 and 2019 growing seasons to devise a new LAI algorithm for crops in the Bekaa Valley. The semi-empirical relationships between the measured LAI and the vegetation indices for all combined crops and different crop groups are given in Table 6. The different types of fitting functions with the crop group-specific LAI-VI relationships are summarized in Table 6 for the two satellites combined (L30 and S30), and in Table 7 for the two satellites separated.
The R 2 for each index obtained with different fittings ranged between 0.27 and 0.63 for the crops combined, with SeLI index having the lowest R 2 in all fittings (0.27). However, the p-value was <0.0001 in all cases. Results for Pearson's r shows a strong correlation between the measured LAI and EVI2 (r = 0.7), followed by SAVI (r = 0.67) and NDVI (r = 0.6) and SeLI (r = 0.5) with the lowest correlation. Each index was represented as a function of the measured LAI for the different crop groups and the combined crops ( Figure 9). However, EVI2 herein can be identified as a best-suited index for a unified algorithm for the crops combined. , and the region between these two vertical lines represents the medium observed LAI range (2-4 m 2 /m 2 ), the 2-D density plots represent four color levels with a dark red color at high density to dark green at lower densities, in addition to the orange and light green in representing medium densities.

Empirical Relationships between VIs and LAI Observations in Bekaa
Vegetation indices are shown to correlate well with LAI. The studied VIs were evaluated with the multi-crop observed LAI dataset over both 2018 and 2019 growing seasons to devise a new LAI algorithm for crops in the Bekaa Valley. The semi-empirical relationships between the measured LAI and the vegetation indices for all combined crops and different crop groups are given in Table 6. The different types of fitting functions with the crop group-specific LAI-VI relationships are summarized in Table 6 for the two satellites combined (L30 and S30), and in Table 7 for the two satellites separated.
The R 2 for each index obtained with different fittings ranged between 0.27 and 0.63 for the crops combined, with SeLI index having the lowest R 2 in all fittings (0.27). However, the p-value was <0.0001 in all cases. Results for Pearson's r shows a strong correlation between the measured LAI and EVI2 (r = 0.7), followed by SAVI (r = 0.67) and NDVI (r = 0.6) and SeLI (r = 0.5) with the lowest correlation. Each index was represented as a function of the measured LAI for the different crop groups and the combined crops ( Figure 9). However, EVI2 herein can be identified as a best-suited index for a unified algorithm for the crops combined.     To evaluate the significance of combining both Landsat-8 and Sentinel-2A in determining the LAI-EVI2 for each crop group, we fitted modeled LAI from L30, from S30, and from both satellites combined using the regression coefficients from Tables 6 and 7 and then extracted the modeled LAI on the observed days. The summary of comparisons between the modeled and observed LAI is shown in Table 8. Our results show that combining both satellites allowed us to achieve lower RMSE and MAE for the crops combined and for each crop group. The use of L30 separately produced higher RMSE and MAE values than those obtained from L30-S30 combined, whereas the use of S30 yielded a lower RMSE and MAE values than L30. Overall, these results demonstrate the benefits of the combined use of Landsat-8 and Sentinel-2A satellites comparing to the single-satellite usage in defining the LAI-EVI2 relationship. Table 8. Statistical comparison between modeled and observed LAI separated by satellite type: L30, S30, and L30 and S30 combined, the modeled LAI is obtained using the LAI-EVI2 empirical relationships displayed in Table 6.

Crop Height
In addition to LAI, crop height is often used to monitor crop development [75][76][77][78][79] and is an essential variable in many land-surface modeling systems-governing surface turbulence and transport of wind and radiation through the canopy to the soil surface [80][81][82][83]. LAI vs. height relationships are unique to specific crop types due to the vast diversity of canopy structure among crop types [84][85][86][87]. In this study, a linear regression analysis was conducted to obtain the relationship between the crop height, the measured LAI, and the best performing LAI models for potato, cannabis, and tobacco. Figures 10-12 indicate that an increase in leaf area is associated with an increased height due to leaf extension as the crop develops. Our analysis shows that the relation between the ground measured LAI and the measured height of the potato crop is essentially linear with an R 2 value of 0.82, greater than that obtained when regressing the crop height to LAI obtained from Model 7 (EVI2, row crops; R 2 = 0.48) (Figure 10). A similar crop height-LAI linear relationship was obtained for the wheat crop against measured LAI with a higher R 2 (0.79) compared to Model 9 LAI (EVI2, wheat)-height relation (R 2 = 0.61) (Figure 11). Furthermore, a good relationship exists between the measured height and the measured LAI for the cannabis and tobacco crops with R 2 : 0.489 and 0.495, and between the modeled LAI (model 10 derived from SeLI), with R 2 : 0.384 and 0.494 for both cannabis and tobacco, respectively ( Figure 12).     To evaluate the importance of combining both L30 and S30 in defining the LAI-h relationship, we applied the L30, S30, and L30 and S30 regression equations obtained between the modeled height and observed LAI for each crop group at each observation day. Statistical comparison between separated and combined satellite datasets is shown in Table 9. This relationship is better defined when both L30 and S30 are combined, where slightly lower RMSE and MAE values were obtained when the LAI based height relationship was determined using both L30 and S30, than those obtained with S30 or L30, separately.

Potato Dry and Fresh Weight
Potato crop biomass is an essential indicator of crop development, and the availability of this information during the growing season can help in decision-making processes for optimal crop production [88,89]. Reflectance-based VI response to canopy variables such as LAI largely determines potato growth, and VIs are widely used in biomass prediction models [11,12]. Here, we regressed the To evaluate the importance of combining both L30 and S30 in defining the LAI-h relationship, we applied the L30, S30, and L30 and S30 regression equations obtained between the modeled height and observed LAI for each crop group at each observation day. Statistical comparison between separated and combined satellite datasets is shown in Table 9. This relationship is better defined when both L30 and S30 are combined, where slightly lower RMSE and MAE values were obtained when the LAI based height relationship was determined using both L30 and S30, than those obtained with S30 or L30, separately.

Potato Dry and Fresh Weight
Potato crop biomass is an essential indicator of crop development, and the availability of this information during the growing season can help in decision-making processes for optimal crop production [88,89]. Reflectance-based VI response to canopy variables such as LAI largely determines potato growth, and VIs are widely used in biomass prediction models [11,12]. Here, we regressed the above-ground dry and fresh weight of the potato crop studied in 2019 against each of the measured LAI, LAI obtained from Model 7 (EVI2, row crops), and the measured height. The relationship between above-ground fresh and dry weights as a function of both the measured LAI and modeled LAI shows that an increase in the leaf area index is associated with an increase in biomass ( Figure 13). No saturation in the measurement of biomass from LAI was noticed. Fresh and dry above-ground biomass of potato could be derived from the modeled LAI using a linear relation was with good accuracy (R 2 : 0.7 for fresh biomass) and with moderate accuracy (R 2 : 0.4 for dry biomass).
On the other hand, the regression of above-ground dry and fresh biomass (AGDB and AGFB, respectively) of potato against crop height showed that an increase in the above-ground fresh and dry weight of the potato crop is associated with an increase in the crop height ( Figure 13). The ability of the linear regression to estimate the fresh and dry above-ground biomass of potato from crop height was with a good accuracy level (R 2 : app. 0.6), with lower RMSE and MAE values obtained as a result of combining both L30 and S30 compared to the single-use of L30 in defining the AGFB and LAI relationship (Table 10). Table 10. Statistical comparison between modeled and observed above-ground dry biomass (AGDB) and above-ground fresh biomass (AGFB) for potato separated by satellite type: L30, S30, and L30 and S30 combined, modeled fresh and dry above-ground biomass is obtained from equations displayed in Figure 13 as a function of observed LAI.

Seasonal Co-Variability in Observed and Modeled Biophysical Variables
This analysis shows that changes in both observed and modeled parameters (LAI, crop height, and biomass) are proportional to changes in EVI2 (Figures 14 and 15). Continuous capture of the biophysical variables' time-series was obtained when combining both L30 and S30 as compared to single sensor use. Modeled LAI derived from the LAI-VI models established in Section 3.2 was used to derive the height and biomass for potatoes on the LAI-height and LAI-biomass equations generated in Section 3.3. We further investigated the seasonal co-variability in the biophysical variables for wheat and potato crops sampled in the Bekaa Valley during the 2019 growing season (Figures 14  and 15). Fresh and dry biomass weight was also evaluated for potato. Measurements of potato crop development were collected regularly during the 2019 growing season, whereas wheat biophysical parameter measurements were conducted at the end of the growing season.
The phenological differences between both potato and wheat were captured for three potato and wheat farms. Typically, the potato had a higher peak of both modeled and measured LAI during the middle of the growing season, where the average maximum measured LAI for all studied potato fields was found to be 6-8 m 2 /m 2 in the mid-season stage (DOY 180). The average minimum measured potato LAI was about 2-3 m 2 /m 2 during the beginning of the growing season (DOY 140). Similarly, modeled LAI evolution was observed for the three potato farms. On the other hand, the measured wheat LAI ranged between 2 and 3 m 2 /m 2 . Lower wheat LAI values were observed since wheat was monitored during the end of the season, where LAI tends to decrease.  single sensor use. Modeled LAI derived from the LAI-VI models established in Section 3.2. was used to derive the height and biomass for potatoes on the LAI-height and LAI-biomass equations generated in Section 3.3. We further investigated the seasonal co-variability in the biophysical variables for wheat and potato crops sampled in the Bekaa Valley during the 2019 growing season (Figures 14 and 15). Fresh and dry biomass weight was also evaluated for potato. Measurements of potato crop development were collected regularly during the 2019 growing season, whereas wheat biophysical parameter measurements were conducted at the end of the growing season.  The phenological differences between both potato and wheat were captured for three potato and wheat farms. Typically, the potato had a higher peak of both modeled and measured LAI during the middle of the growing season, where the average maximum measured LAI for all studied potato fields was found to be 6-8 m 2 /m 2 in the mid-season stage (DOY 180). The average minimum measured potato LAI was about 2-3 m 2 /m 2 during the beginning of the growing season (DOY 140). Similarly, modeled LAI evolution was observed for the three potato farms. On the other hand, the measured wheat LAI ranged between 2 and 3 m 2 /m 2 . Lower wheat LAI values were observed since wheat was monitored during the end of the season, where LAI tends to decrease.
The measured crop height increases as the crop develops. The average minimum measured potato crop height was about 0.4 m during the beginning of the growing season (DOY 140). The maximum average crop height for all studied potato fields was found to be around 0.6 m in the midseason stage (DOY 180). The measured wheat crop height was found to be approximately 1 m during the end of the season (DOY 135-140). Similar modeled crop height evolution was observed for the three wheat farms.
Similarly, the above-ground fresh and dry biomass increased with crop growth. The aboveground biomass was the lowest at the beginning of the growing season, and it increases to reach its Similarly, the above-ground fresh and dry biomass increased with crop growth. The above-ground biomass was the lowest at the beginning of the growing season, and it increases to reach its maximum during the middle of the growing season with a value of 2 kg/m 2 and 0.2 kg/m 2 , for the above-ground fresh and dry biomass, respectively. Different temporal LAI, crop height, and biomass dynamics were noticed in both potato and wheat crops (Figures 14 and 15). Both crop type and management practices affect the dissimilarity of the temporal biophysical parameters of these crops. Furthermore, the same crop types, but different fields, diverged in their biophysical temporal profiles. These deviations can be attributed to the variation of farming conditions, agronomic, and management practices among the farms with the same crops. Figure 16 shows the time development of the LAI, crop height, and above-ground fresh and dry biomass of a potato crop surveyed during the 2018 season. The modeled LAI maps are derived from the exponential LAI-EVI2 fit established in Section 3.2. These LAI maps were used to derive the height and biomass for potatoes on the LAI-height and LAI-biomass equations generated in Section 3.3. The spatial variability of the canopy conditions throughout the growing season is considerable, reflecting differences in soil conditions, planting dates, and cultivation practices. Note that fields 10 and 11 have depressed growth patterns coinciding with the L30 satellite overpass (DOY: 167). Furthermore, fields 9, 10, 11 and half of 8 reached senescence earlier than the other fields. The fact that fields 10 and 11 lag behind other potato-fields in terms of crop development is probably due to soil moisture stress that could be an insufficient amount of water available for irrigation or that irrigation was delayed. However, moisture stress was not high enough to completely inhibit crop growth.

Spatial Biophysical Parameter Patterns
Other factors yielding the variation in the magnitude of biophysical parameters estimation include the crop's susceptibility to unfavorable growth conditions, crop variety, differences in planting and harvest dates, crop management, weed control, planting density, and irrigation management. Overall, HLS products present good spectral and temporal resolution allowing for visual inspection of error estimates in time series maps of LAI and other biophysical variables.

Uncertainties in -LAI Assessment
Uncertainties associated with LAI assessment can be grouped into two categories: the uncertainty of remote sensing methods and the errors in the LAI observation techniques.
LAI can be retrieved using a physical radiative transfer model or a simple empirical model based on a vegetation index. Any physical model is a simplification of the complex world. There are uncertainties associated with the physical model itself, and in the surface reflectance data used as input. The uncertainty of atmospheric correction of the satellite images is a major factor affecting the accuracy of LAI estimation by remote sensing [58]. Atmospheric correction also affects LAI-VI relationships. The use of vegetation indices to generate LAI products is subject to several uncertainties. First, studies have shown that due to VI saturation effects, LAI values that are greater than 4 m 2 /m 2 are beyond the prediction power of the studied vegetation indices [5,24,60]. Thus, issues with VI saturation over dense canopies remain an issue to be addressed [58,90]. Second, different vegetation indices respond differently to the crop leaf optical properties defined by the leaf chlorophyll, dry mass, and water content [67]. Third, the leaf inclination, solar angle, and sensor bandwidth and calibration affect the vegetation spectral reflectance [91]. Therefore, LAI-VI relationships usually are only good for the local area where model was developed, while physical models could, in theory, be used more generally.
Regarding errors in the LAI measurement method, the underestimation of LAI in the high range of measured LAI can be attributed to the SS1 equipment limitations under high vegetation cover (LAI range > 4 m 2 /m 2 ). LAI that is calculated using light transmission theory, such as the SS1 method, is underestimated because SS1 does not correct for canopy clumping effects [92]. Generally, when the canopy is clumped, the leaves hide each other, leaving a gap in between, and thus allow for more light to reach the ground as compared to randomly distributed leaves. The underestimation of LAI with SS1 was also observed in previous studies due to the lack of clumping correction [93,94]. Another type of error that falls in the second uncertainty category is that the leaf angle distribution (ELADP) is a fixed constant that is predefined by the user according to the crop. However, the chosen ELADP values may not be accurate throughout the growing season since the canopy leaves change throughout the entire season [95]. The assumption of a constant leaf absorption constant of 0.85 adds uncertainty to the LAI measurements. Moreover, row spacing, crop height, and time of measurement can affect the LAI values.
Scaling up field observations to match satellite pixel footprint is another challenge. Field LAI is a point measurement. Even through many points can be measured within a satellite pixel, the summary of the point measurements is still an approximation of the area. In addition, satellite pixels have geolocation errors and their footprints may vary from date to date, which is always a challenge for remote sensing data product validation [96].

S2 Model Results in Comparison to Previous Studies
We also assessed the performance of the newly developed SeLI-based LAI algorithm. Pasqualotto et al. [30] showed that LAI derived SeLI could be used as a generic model for LAI green of different crop types. Our results agree with the conclusions obtained by Pasqualotto et al. [30] and Sharma, et al. [97], who pointed out that the S2 red-edge bands are critical for the estimation of biophysical variables, mainly the LAI. Following the HLS VI models, Model 10 (SeLI, developed for overall agricultural crops) ranked next in terms of low RMSE, MAE, MBE, RRMSE, and MAPE values for all crops sampled in the study area. However, the S2 LAI algorithm significantly underestimates high LAI compared to other LAI-HLS derived models.
Our assessment shows that LAI derived from SNAP performed best for grain crops. The results demonstrate the potential of using SNAP to derive more uniform LAI estimates for both wheat and barley (possibly due to the good performance of the radiative transfer model for homogenous conditions). A clear underestimation of the SNAP-LAI product was noticed for herbs, potato, vegetables, and tobacco ( Figure 8). Pasqualotto et al. [30] showed that LAI estimated from the Sentinel-2 product using the SNAP software underestimated LAI of potato, artichoke, squash, alfalfa, lettuce, wheat, pumpkin, and others. Compared to HLS VI and S2 SeLI models, SNAP provided less accurate estimates of LAI (Table 5 and Figure 8). These results are not in agreement with Kganyago et al. [98] who showed that SNAP-derived LAI values were overestimated for maize and sunflowers during the 2019 peak growing time in South Africa. This shows that the LAI estimation using SNAP software is highly location and crop type dependent.
S2 and HLS models outperforming the SNAP model confirm the suitability of S2 models applications to the small geographical scale (field level) of our study. The authors of the SNAP (personal communication) mentioned that the biophysical processor lacks the ancillary information from which some vegetation architecture features could be derived. Hence, the absence of reliable and regularly updated land cover maps at a spatial resolution, like that of S2, prevents tailoring of specific algorithms for each of the land cover classes. Consequently, the algorithm may provide reasonably good estimates over all the cases, but certainly poorer performances as compared to algorithms specific to a given surface type [28]. Thus, the poor performance of SNAP-derived LAI suggests that it is not suitable for precision agriculture and may be useful for large-area applications.

HLS VI Models Validation Challenges
The major challenge encountered while validating existing vegetation index relationships applied to HLS SR product is the uniqueness and diversity of the site-specific LAI-VI relationships. Site-specific properties such as the crop type [5], biochemical properties of the leaf canopy [99], biophysical properties of the leaf canopy [100], as well as sensor characteristics and the atmospheric scattering and absorption [101], have yielded unique local LAI-VI relationships that are ideally applicable only under unique conditions. This work shows that the existing LAI-VI models, initially developed for particular crops from different sensors, can be applied to some crops sampled in this research irrespective of their limitation to "one place, one time, one equation." As shown in Results, the LAI models originally developed for wheat (RMSE: 1.27), maize (RMSE: 1.34), and row crops (RMSE: 1.38) worked well for the combined crops of the study area. However, these relationships cannot replace the locally established relationships in this work, as the local relationships account for the site-specific conditions affecting the LAI-VI relationships.

Empirical Relationships Establishment between VIs and LAI Observations Challenges
The major challenge of establishing the empirical relationship between the observed LAI and vegetation indices is the difficulty of finding one index with a sufficiently general character that can be used to estimate LAI for a wide variety of crop types. In this work, the established indices do not present this general character, as each index performed differently in different crop types. The most robust relation between ground measured LAI and vegetation indices appears in SeLI for cannabis, SAVI and EVI2 for potato, NDVI for tobacco, and EVI2 for vegetables. However, EVI2 outperformed other vegetation indices for all crops combined, and therefore it can be identified as an index that is best suited for a unified algorithm for crops in semi-arid irrigated regions with heterogeneous landscapes.
As more remote sensing and field measurements become available, machine learning approaches may be used to investigate the complex relationships between LAI and remote sensing spectral features. Surface reflectance provides information about land covers and structures, which could be used to build the relationship directly [102].

Conclusions and Recommendations
This study reconfirms the limited ability of Landsat as an isolated sensor to capture the full time-series at the farm level. L30 missed biophysical observations at critical growth stages due to the smaller number of images available. Few L30 images were used due to its long revisit time and cloud coverage. Our results also demonstrate the importance of higher observation frequency achieved with the combination of Landsat-8 and Sentinel-2A satellites compared to the single use of L30 in both crop monitoring and modeling efforts, where stronger LAI-based crop biophysical parameters relationships were achieved with L30 and S30 combined, as compared to L30 separately. LAI and other crop growth variables are appropriately monitored and modeled in this work, especially during fast-developing phenomena. Multi-satellite data normalization is crucial, and this was addressed through crop biophysical variables time-series generated for monitoring wheat and potato crops in the Bekaa Valley. Time series of crop growth parameters are essential for detecting changes in crop cover, and the combination of data from different ongoing satellite missions provides this opportunity. Thus, the recent availability of high temporal frequency high-resolution images (HLS product) provides the ability to conduct and utilize more site-based biophysical variables measurements for crops other than the two studied strategic crops (wheat and potato) by providing more time-matched surface reflectance products.
Using the combined Landsat 8 and Sentinel-2 data at a 30-m spatial resolution, we focused in this work on the estimation LAI at the field scale in a semi-arid irrigated heterogeneous landscape. Unlike most existing approaches utilizing vegetation indices derived from non-normalized satellite data to estimate LAI, we derived the LAI from three methods, namely: existing vegetation index (VI) relationships applied to Harmonized Landsat-8 and Sentinel-2 (HLS) surface reflectance produced by NASA, the SNAP biophysical model, and L2A THEIA product from Sentinel-2. The EVI2 index surpassed any other index in predicting LAI, and it is the best for modeling LAI of various crops. This study also shows that LAI derived from the artificial neural network through ESA's SNAP biophysical processor is underestimated but yielded a good accuracy level for wheat and barley LAI estimation. Hence, the poor performance of SNAP dictates further improvement to support precision agriculture [98]. Future studies should investigate the discrepancy in the disagreement of LAI-SNAP accuracy in the cited research work. While the LAI derived from the SeLI algorithm produces an acceptable accuracy level compared to HLS-EVI2, a significant underestimation at high LAI was observed. The normalizing satellite data in crop monitoring and modeling is crucial to eliminate the bias due to sensor and bandwidth variability and atmospheric correction techniques. We encourage the scientific community to validate new crop and site-specific LAI-VI, LAI-height, and LAI-above-ground biomass empirical relationships developed in this work.
We recommend exploring improved hemispherical LAI estimation techniques that can provide information on the canopy clumping and leaf angle inclination, measure the LAI at different heights, and differentiate between leaves and woody parts through the integration of infrared techniques [103,104]. The use of other vegetation indices, such as the normalized difference water index (NDWI), derived from near-and shortwave-infrared reflectance, may be less susceptible to saturation problems [105]. Ingestion of multiple spectral bands within a machine learning approach may be useful for large-area applications [102]. Future work may also focus on the effects of soil moisture on the retrieval of LAI. More studies on the use of other vegetation indices that correct for soil background, such as the Weighted Difference Vegetation Index (WDVI) [106], could be explored (here we found that EVI2, which corrects for both atmosphere and soil background, gave good results). Additionally, future studies should consider quantifying the effect of various atmospheric correction approaches such as Sen2Cor, FORCE, iCor, and MAJA on biophysical parameters retrieval, since atmospheric correction performance varies by different spectral bands of Sentinel-2 [98].

Acknowledgments:
We thank Naji El-Beyrouthy for helping with field data collection.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
There are two major classes of LAI measurement methods: (1) direct measurement and (2) indirect measurement. The direct method involves destructively harvesting the canopy and measuring the area of the collected leaves samples using instruments such as LI-3000C portable leaf area meter (LI-COR, Lincoln, NE, USA). This method has relatively high precision, and is often used to validate the indirect measurement methods; however, it is very labor-intensive, tedious, and is not suitable for large area measurement [107,108]. On the other hand, the indirect methods are based on light interception, hemispherical photography, or remote sensing [15]. Both light interception and hemispherical techniques are referred to as optical approaches, which involve variable measurements such as gap fraction and light transmission. These variables are measured by commercial instruments and are related to LAI based on PAR inversion and Beer-Lambert extinction law theories for both gap fraction and light transmission, respectively. A more detailed explanation of these theories is provided in Appendix A, along with an evaluation of different optical instruments as presented in literature.
Several studies have used different optical instruments to either assess these instruments against direct measurements or to validate the remotely sensed LAI [15,95,109]. Although indirect measurement of LAI using optical instruments is faster than direct measurement, it suffers from several limitations. First, indirect measurements do not measure the LAI but rather the light intercepted by the plant canopy that is then related to LAI. In other words, the optical instruments do not distinguish the leaves from the flowers or branches. Thus, they tend to either underestimate or overestimate the direct LAI depending on the method used and canopy studied.
Moreover, each method has sources of error. For example, the digital hemispherical photography (DHP) includes some subjectivity regarding factors such as variable clouds, shadows, and brightness, which makes it challenging to choose a threshold that allows the distinction of the canopy from non-canopy. For this reason, hemispherical photographs are recommended to be captured under uniformly diffuse or overcast conditions, while the leaf angle distribution can be a source of error if it is misestimated in instruments such as SS1 and AccuPAR. Finally, row spacing, crop height, time of measurement, and placement of the meter can affect the observed LAI values.
The indirect measurement of LAI with the optical instruments is based on two major theoretical approaches, namely: the gap fraction or PAR inversion methods.
With the gap fraction approach, LAI is determined by solving for L as per the below Equations (A1) and (A2): where P (θ) is the gap fraction, θ is the zenith angle of view, α is the leaf angle, G (θ, α) is the fraction of foliage projected on the plane normal to the zenith direction, and Ω (θ) is the clumping coefficient when foliage is randomly distributed within the canopy, Ω (θ) = 1, but as the foliage becomes more clumped, Ω (θ) < 1, and L is the leaf area index. Examples of instruments that use the above theory include the LAI-2200 plant canopy analyzer, which compares above and below canopy light levels detected by five conical rings with zenith view angles ranging from nadir to 75 • . Another example is digital hemispherical photography (DHP). DHP measures LAI by the inversion of the Poisson model (Equation (A2)), as described by Thimonier, et al. [110]. DHP uses a fish-eye lens inside the camera that images the canopy from beneath in a hemisphere. The collected photographs are then processed by the user using Can-Eye or other software. By process of thresholding, distinguishing pixels that are occupied by leaves from those that are occupied by the sky, the user can determine the LAI values.
The second approach is based on Beer-Lambert extinction law, and it can be explained as follows: In a case of light absorption by the canopy, according to Beer's law, a relationship between the incident (I 0 ), transmitted light (I) and leaf area index (L) is given by Equation (A3) below: where K is an extinction coefficient which depends on the direction of beam and leaf angle distribution. K is calculated, according to Campbell [111] as per Equation (A4): where x is the ellipsoidal leaf angle distribution (ELADP), and θ is the zenith angle of the direct beam.
Typical ELADP values used in this study are summarized in Table A1 from Campbell and Van Evert [112]. A default value of 1 was used for random spherical distribution crops that are not mentioned in Table A1 below. Zenith angels are calculated from longitude, latitude, and local timing. Examples of instruments that use the above theory include SS1 SunScan Plant Canopy Analysis System. This instrument consists of 64 PAR sensors with an accuracy of ±10%, embedded in a 1 m long portable probe. Combined with a calibrated BF3 Sunshine Sensor, it measures the incident and transmitted PAR in the canopy to provide an estimate of LAI. A brief theoretical background of how the SunScan computes its reading of LAI can be found in the user's manual. Another instrument is the Accu-PAR model Lp-80 PAR LAI Ceptometer (Decagon Devices, Inc. Pullman, WA, USA). The LAI measured by Accu-PAR and other PAR inversion instruments is dependent on several factors such as the solar zenith angle which is calculated using the time of day and the geographic location, the path length of the beam radiation, and the leaf angle distribution of the canopy as described by Beer's law.
Evaluation of different LAI instruments as presented in Literature: Casa et al. [113] evaluated different commercial LAI instruments and compared them to direct LAI measurements of different crops. These instruments are the LAI-2000, the Sunscan Ceptometer, and digital hemispherical photography. Casa et al. [113] concluded that the digital hemispherical photograph provided the best estimates of LAI for alfalfa (RMSE = 0.33 m 2 /m 2 ), broad bean (RMSE = 0.3 m 2 /m 2 ), emmer (RMSE = 0.31 m 2 /m 2 ), maize (RMSE = 0.58 m 2 /m 2 ), and wheat (RMSE = 0.68 m 2 /m 2 ). Ariza-Carricondo et al. [15] assessed the LAI-2200 plant canopy analyzer, the SS1 SunScan Canopy Analysis System, and digital hemispherical photography (DHP) for four different canopy types, including maize. The results have shown that LAI-2200 overestimates LAI in a very low range of the directly measured LAI. The DHP only slightly overestimated the direct method, whereas the SS1 underestimated the direct measurements by about 25%; this is consistent with Gower et al. [94], who showed that SunScan underestimated LAI compared with the destructive measurement. However, the SunScan probe has provided reasonable estimates of leaf area index for corn (p = 0.0713) as reported by Wilhelm et al. [95] and for rice cultivars (p < 0.001) when restricting the range of LAI to <4 m 2 /m 2 as reported by Sone et al. [110].
Reyes-González et al. [3] measured the leaf area index of cornfields with the Accu-PAR model Lp-80 PAR LAI Ceptometer (Decagon Devices, Inc. Pullman, WA, USA). In their methodological approach, the LAI measurement was done during the time the satellite overpasses to assess the LAI estimated using the remote sensing-based METRIC model. The results of the comparison of the estimated LAI to the LAI obtained from AccuPAR have shown that LAI values measured in situ were higher than the LAI values estimated using the remote sensing-based METRIC model by about 12%, with a good linear correlation (MBE of 0.61, and RMSE of 0.59). Pasqualotto et al. [30] and Anderson et al. [105] used LAI-2200 plant canopy. Pasqualotto et al. [30] compared LAI obtained from Sentinel-2 LAI Index SeLI ([B8a − B5]/[B8a + B5]) with the in situ measurements and found a good linear correlation between them (R 2 of 0.732, and RMSE of 0.69 m 2 /m 2 ). Similar results were achieved by Anderson et al. [105], where the accuracy of LAI retrieval from optimized soil adjusted vegetation index (OSAVI) and normalized difference water Index (NDWI) compared to the measured LAI was 0.6 as measured by the root-mean-square-deviation (RMSD).