Using the Negative Soil Adjustment Factor of Soil Adjusted Vegetation Index (SAVI) to Resist Saturation Effects and Estimate Leaf Area Index (LAI) in Dense Vegetation Areas

Saturation effects limit the application of vegetation indices (VIs) in dense vegetation areas. The possibility to mitigate them by adopting a negative soil adjustment factor X is addressed. Two leaf area index (LAI) data sets are analyzed using the Google Earth Engine (GEE) for validation. The first one is derived from observations of MODerate resolution Imaging Spectroradiometer (MODIS) from 16 April 2013, to 21 October 2020, in the Apiacás area. Its corresponding VIs are calculated from a combination of Sentinel-2 and Landsat-8 surface reflectance products. The second one is a global LAI dataset with VIs calculated from Landsat-5 surface reflectance products. A linear regression model is applied to both datasets to evaluate four VIs that are commonly used to estimate LAI: normalized difference vegetation index (NDVI), soil adjusted vegetation index (SAVI), transformed SAVI (TSAVI), and enhanced vegetation index (EVI). The optimal soil adjustment factor of SAVI for LAI estimation is determined using an exhaustive search. The Dickey-Fuller test indicates that the time series of LAI data are stable with a confidence level of 99%. The linear regression results stress significant saturation effects in all VIs. Finally, the exhaustive searching results show that a negative soil adjustment factor of SAVI can mitigate the SAVIs’ saturation in the Apiacás area (i.e., X = −0.148 for mean LAI = 5.35), and more generally in areas with large LAI values (e.g., X = −0.183 for mean LAI = 6.72). Our study further confirms that the lower boundary of the soil adjustment factor can be negative and that using a negative soil adjustment factor improves the computation of time series of LAI.


Introduction
As major land cover of the planet, forests have become a key priority in studies of the biodiversity and the carbon cycle of terrestrial ecosystems [1,2]. It is essential to record forest dynamics to understand the terrestrial carbon cycle better and improve forest management practices [3,4]. Leaf area index (LAI), defined as the one-sided green leaf area per unit ground area in broadleaf canopies and one-half the total needle surface area per unit ground area in coniferous canopies, is an essential indicator for describing the canopy structure of forest ecosystems, estimating the primary productivity of the stand, and evaluating forest condition over large areas [5,6].
Our study uses two data sets: MODIS 4-Day global 500 m LAI product (MCD15A3H V6 level 4, provided by the National Aeronautics and Space Administration Land Processes Distributed Active Archive Center at the U.S. Geological Survey Earth Resources Observation and Science Center which is located outside Sioux Falls, South Dakota, USA) from time-series observation data and the global LAI product global leaf area index from field measurements (GLAIFM) from field measurements. MODIS LAI data were used to study a region, Apiacá, for which VIs were calculated from Sentinel-2 and Landsat-8 surface reflectance. For the global field-measured LAI dataset, the VIs were calculated with Landsat-5 surface reflectance products. global 500 m LAI products (MCD15A3H V6 level 4). They were used to validate LAI estimation from VIs calculated from surface multispectral reflectance datasets computed with Sentinel-2 images from 4 February 2019, to 10 November 2020, and Landsat-8 Tier 1 images from 16 April 2013, to 21 October 2020. The use of multiple sources of data helped us to increase the sampling frequency of the time series. Because of the similarity of satellite resolution of Sentinel-2 (10 m) and Landsat-8 (30 m), the VIs differences caused by the two sensors were neglected. First, the multispectral images were processed. All available Sentinel-2 and Landsat-8 Tier 1 surface reflectance data in GEE before 10 November 2020, were selected and then filtered by two rules: selecting images that cover the study area and an average cloud cover lower than 20%. Then, each remaining image was clipped to ensure that only the study area's pixels were retained. After that, two boundary masks were used in a preprocessing stage. The first was the cloud mask, and the second was the vegetation mask. GEE provides cloud masks of Sentinel-2 and Landsat-8. Pixels in the cloud mask were removed in the subsequent data processing. The vegetation mask of Sentinel-2 was created through the classification of the NDVI and simple ratio NIR/Green Ratio Vegetation Index (GRVI) [29] images using thresholds indicated in a technical report of Sentinel-2 [30]: pixels with (NDVI > 0.40) or (0.36 ≤ NDVI ≤ 0.40 and GRVI > 2.50) were marked as vegetation pixels. The same rule was implemented to create a Landsat-8 vegetation mask. VIs and reflectance values were calculated pixelwise, and the final VIs and reflectance of each satellite The multispectral satellite datasets and MODIS LAI datasets used in this study were provided and processed by the Google Earth Engine [26][27][28] (GEE, https://earthengine. google.com/, accessed on 11 March 2021). LAI data were from the MODIS 4-Day global 500 m LAI products (MCD15A3H V6 level 4). They were used to validate LAI estimation from VIs calculated from surface multispectral reflectance datasets computed with Sentinel-2 images from 4 February 2019, to 10 November 2020, and Landsat-8 Tier 1 images from 16 April 2013, to 21 October 2020. The use of multiple sources of data helped us to increase the sampling frequency of the time series. Because of the similarity of satellite resolution of Sentinel-2 (10 m) and Landsat-8 (30 m), the VIs differences caused by the two sensors were neglected.
First, the multispectral images were processed. All available Sentinel-2 and Landsat-8 Tier 1 surface reflectance data in GEE before 10 November 2020, were selected and then filtered by two rules: selecting images that cover the study area and an average cloud cover lower than 20%. Then, each remaining image was clipped to ensure that only the study area's pixels were retained. After that, two boundary masks were used in a preprocessing stage. The first was the cloud mask, and the second was the vegetation mask. GEE provides cloud masks of Sentinel-2 and Landsat-8. Pixels in the cloud mask were removed in the subsequent data processing. The vegetation mask of Sentinel-2 was created through the classification of the NDVI and simple ratio NIR/Green Ratio Vegetation Index (GRVI) [29] images using thresholds indicated in a technical report of Sentinel-2 [30]: pixels with (NDVI > 0.40) or (0.36 ≤ NDVI ≤ 0.40 and GRVI > 2.50) were marked as vegetation pixels. The same rule was implemented to create a Landsat-8 vegetation mask. VIs and reflectance values were calculated pixelwise, and the final VIs and reflectance of each satellite image were derived by aggregating VIs and reflectance calculated per pixel. In the next step, LAI was extracted for each corresponding satellite image.
In this work, we used MODIS leaf area index 4-day global 500 m (MCD15A3H V6, https://doi.org/10.5067/MODIS/MCD12Q1.006, accessed on 11 March 2021) [31]. LAI data were searched for each surface reflectance image within eight days of the acquisition time of the surface reflectance image. If no LAI product was available, LAI data within 16-day intervals was searched instead. If no LAI data was available within 16 days, this particular surface reflectance image was ignored. After getting the LAI data in the time range of satellite images, we extracted pixels within both the study area boundary and reflectance image boundary, and a mean LAI value was derived by aggregating all pixels' LAI inside these two boundaries. Finally, for each multispectral image, we had a mean value of different VIs, a mean value of reflectance in various bands, and a mean value of LAI.
Depending on the type and location of sensors, the study area could correspond to several remote sensing images on the same day. In that case, mean values of VIs, reflectance, and LAI were calculated using the images acquired on that day. The reason for considering the entire study area rather than pixels as one study object was to smooth out inconsistencies caused by differences in the resolution of the remotely sensed data.

Global Field Measured LAI Data
We used global LAI data from the GLAIFM dataset [32]. This dataset is compiled from around 1000 published estimates of LAI covering the period from 1932 to 2000. These historical LAI data include natural and seminatural (managed) ecosystems and some cultivated vegetation. It contains a wide range of LAI values of 15 biome/land cover classes, from 0.46-2.16 for deserts to 4.4-13.04 for tree plantations. We downloaded it from the website (https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=584, accessed on 11 March 2021).
We extracted the Landsat-5 Tier 1 surface reflectance product for each LAI pixel using the LAI measurement date and geographical information (latitude and longitude). Then, we calculated VIs from these reflectance values. The specific steps are as follows: (1) Data cleaning: data with null values for latitude, longitude, or date of measured LAI were removed. Also, data with LAI as a range value instead of a specific value were removed; data with LAI measurements spanning over one month were removed. (2) Screening of data corresponding to dense vegetation: data with LAI values less than 4.0 were removed. (3) Screening of time range: data with a LAI measurement date range outside the Landsat-5 coverage time range were removed. (4) Removing duplicate data: some data had the same measurement time and geographical location with different LAI values. On closer examination, we found that some of these data were from measurements of different biomes. Therefore, these data were also removed as we could not establish a one-to-one relationship between LAI and VIs. (5) Matching LAI with satellite reflectance image: because the temporal accuracy of all measured LAI data is only one month, we searched all surface reflectance data for the month of Landsat-5 Tier 1 data and averaged them to obtain the final surface reflectance data. Some reflectance data could not be retrieved due to cloud cover, so their corresponding LAI data were removed. Table 1 shows the LAI data and corresponding station information that remained after the above steps. These data were used for evaluating the LAI estimation accuracy of VIs. The LAI value ranges from 4.06 to 10.59, and the sites are mainly located in Canada, Japan, and the USA.
where R and N are the atmospherically corrected surface reflectance in the red and NIR bands.
Conversely to NDVI, SAVI and TSAVI have a soil adjustment factor designed to mitigate the soil noise by considering multiple soil scattering. Here, this factor is referred to as "X" with the VIs name as a suffix. SAVI is defined as: where X SAVI is the soil adjustment factor of SAVI. A value equal to 0.5 is recommended in its original paper [11]. TSAVI is defined as: where a and b are the slope and interception of the soil line, respectively. Here these two parameters were set as the constant value of 1.2 and 0.04, respectively, which are considered global soil line parameters [12]. X TSAVI is the soil adjustment factor of TSAVI, recommended to be equal to 0.08 in the original paper [12]. EVI, in addition to a soil adjustment factor, uses atmosphere resistance factors: where B is the atmospherically corrected surface reflectance in the blue band. G is a factor gain equal to 2.5. C 1 , C 2 are the aerosol resistance coefficient (C 1 = 6, C 2 = 7.5), using the blue band reflectance to correct aerosol influences in the red band. X EVI is a soil adjustment factor initially equal to 1 [19]. Since the surface reflectance products were used as multispectral images, we assumed that atmospheric influences had been almost completely removed. Therefore, atmosphere resistance was not considered in this paper.

Time Series Analysis
Time series analysis is a statistical method that processes time series data to forecast, control, and understand features of the data. [33]. Time series of LAI can be separated into the trend, seasonality, and residuals. Trend is the increasing (or decreasing) value in the series data. Seasonality is the repeating short-term cycle in the data series. Residuals are the Sensors 2021, 21, 2115 6 of 14 time series data after the trend, and seasonal components are removed. The decomposition was achieved using classical seasonal decomposition by moving averages [34]. The seasonal multiplication component was selected, and the base-level static analysis was conducted with the rolling mean and Dickey-Fuller test. As the threshold for this analysis, a 99% significance level was employed for the null hypothesis (i.e., no trend). This test explored the trend from a quarter of a year to the next. In the time series analysis, the independent variable was date, and the dependent variable was VIs or LAI. Considering that the field measured LAI data were not continuously observed, only MODIS LAI data were used for time series analysis.

Linear Regression
Linear regression is a method to build the linear relationship between a dependent variable and one or more independent variables [35]. We used a unitary linear regression with VI as the independent variable and LAI as the dependent variable. The linear regression equation is defined as: where c is the slope and d is the interception of the regression line. The linear regress is widely used to estimate LAI from VIs for its simplicity [8,36]. When LAI increases to a certain threshold value and no longer maintains a linear relationship with VIs, the saturation effect is observed.

Exhaustive Search
For discrete problems where no effective solution method is known, it may be necessary to test every possibility sequentially to determine whether it is the solution. Such exhaustive examination of all possibilities is called exhaustive search, direct search, or the "brute force" method [37]. The exhaustive search was used to evaluate the performance of SAVI calculated with exhaustive optimal soil adjustment factor X SAVI . Here the independent variable was X SAVI , and the dependent variable was R 2 , slope, and P of the linear regression equation, respectively. The work of Ren et al. [38] demonstrated that a negative X SAVI value (e.g., −0.2) is acceptable, which allowed us to set the exhaustive interval of X SAVI as [−0.3, 1] and to divide the boundary into 2000 parts at intervals of 0.001. The ability to estimate LAI using SAVI was investigated using a X SAVI sequence over the whole study period. The estimation was achieved using a linear regression model between LAI and SAVI calculated with the X SAVI sequence. Both MODIS and field measured LAI data were used here. Finally, the R 2 and P were used to evaluate the linear regression performance of various SAVI values calculated with positive or negative X SAVI .

Seasonal Decomposition of LAI
Time series analysis of LAI was conducted to ensure the stability of vegetation conditions in the study area. The result of the decomposition analysis of LAI time series can be seen in Figure 2. All LAI data (labeled as observation in the figure) were averaged using a quarter (i.e., winter, spring, summer, autumn). However, due to cloud cover, few quarters had no corresponding LAI. In that case, a mean value between the preceding and the following values of this null value was assigned. All quarterly LAI data were then decomposed into three parts: trend, seasonal, and residual with a quarterly frequency.
The residual LAI stationarity was validated using the rolling mean ( Figure 3) and the Dickey-Fuller test ( Table 2). As shown in Figure 3, the means and variances of LAI tended to be constant with little volatility. According to the Dickey-Fuller Test results (Table 2), the statistic value was much smaller than the critical value at 1%, meaning that the data can be considered stable at the 99% confidence level. The residual LAI stationarity was validated using the rolling mean (Error! Referenc source not found.) and the Dickey-Fuller test (Error! Reference source not found.). A shown in Error! Reference source not found., the means and variances of LAI tended be constant with little volatility. According to the Dickey-Fuller Test results (Error! Refe ence source not found.), the statistic value was much smaller than the critical value at 1% meaning that the data can be considered stable at the 99% confidence level.

Time Series Analysis of LAI and VIs
The variations of VIs and LAI time series are plotted in Figure 4. It allows one to compare the performance of VIs on LAI estimation. The VIs were calculated with the default soil adjustment factor (X SAVI = 0.5, X TSAVI = 0.08, and X EVI = 1). LAI values ranged Sensors 2021, 21, 2115 8 of 14 from 2 to 6.5. NDVI had the highest value over the entire period (i.e., values from 0.7 to 0.9), followed by EVI (i.e., values from 0.4 to 0.7). SAVI and TSAVI shared almost the same range, from 0.2 to 0.6. These high values suggest that the area has high foliage coverage. Besides, Figure 4 shows that VIs can be indicators or not of LAI's trends depending on time. The solid circles highlight a period when LAI has a trend opposite to that of all VIs. The slash circles highlight a period when the trend of LAI is consistent with NDVI while being opposite to other VIs. The dashed circles highlight a period when the trend of LAI is consistent with all VIs.

Time Series Analysis of LAI and VIs
The variations of VIs and LAI time series are plotted in Error! Reference source not found.. It allows one to compare the performance of VIs on LAI estimation. The VIs were calculated with the default soil adjustment factor (XSAVI = 0.5, XTSAVI = 0.08, and XEVI = 1). LAI values ranged from 2 to 6.5. NDVI had the highest value over the entire period (i.e., values from 0.7 to 0.9), followed by EVI (i.e., values from 0.4 to 0.7). SAVI and TSAVI shared almost the same range, from 0.2 to 0.6. These high values suggest that the area has high foliage coverage. Besides, Error! Reference source not found. shows that VIs can be indicators or not of LAI's trends depending on time. The solid circles highlight a period when LAI has a trend opposite to that of all VIs. The slash circles highlight a period when the trend of LAI is consistent with NDVI while being opposite to other VIs. The dashed circles highlight a period when the trend of LAI is consistent with all VIs.

Linear Regression between VIs and LAI
We tested the ability of VIs for LAI estimation. Error! Reference source not found. shows the scatter plot of LAI (y-axis) and VIs (x-axis). The linear regression model was used to compare the VIs. Error! Reference source not found. gives the coefficients of the

Optimal Soil Adjustment Factor X 4.1. Linear Regression between VIs and LAI
We tested the ability of VIs for LAI estimation. Figure 5 shows the scatter plot of LAI (y-axis) and VIs (x-axis). The linear regression model was used to compare the VIs. Table 3 gives the coefficients of the linear regression models. Figure 5 and Table 3 clearly show strong saturation effects for all VIs. The VIs hardly showed an upward trend with the increase of LAI. Besides, all VIs showed a weak correlation relationship with LAI. In general, NDVI showed relatively good performance in both cases, with the highest R 2 (0.1632 in MODIS LAI and 0.4313 in field measured LAI) and the lowest P (0.0173 in MODIS LAI and 0.2860 in field measured LAI), followed by TSAVI. EVI and SAVI performances were unexpected, especially with the MODIS LAI data, with which they were negatively correlated. In all, the use of VIs to estimate LAI in very densely vegetated areas could lead to large errors. linear regression models. Error! Reference source not found. and Error! Reference source not found. clearly show strong saturation effects for all VIs. The VIs hardly showed an upward trend with the increase of LAI. Besides, all VIs showed a weak correlation relationship with LAI. In general, NDVI showed relatively good performance in both cases, with the highest R 2 (0.1632 in MODIS LAI and 0.4313 in field measured LAI) and the lowest P (0.0173 in MODIS LAI and 0.2860 in field measured LAI), followed by TSAVI. EVI and SAVI performances were unexpected, especially with the MODIS LAI data, with which they were negatively correlated. In all, the use of VIs to estimate LAI in very densely vegetated areas could lead to large errors.

LAI Estimation Using SAVI with Negative Soil Adjustment Factor XSAVI
We also conducted a linear regression analysis of SAVI and LAI by varying XSAVI from −0.3 to 1 with a step of 0.001 to validate SAVI's robustness using a negative soil adjustment factor. Both MODIS and field measured LAI data were used for the validation. The results of the linear regression are shown in Error! Reference source not found..

LAI Estimation Using SAVI with Negative Soil Adjustment Factor X SAVI
We also conducted a linear regression analysis of SAVI and LAI by varying X SAVI from −0.3 to 1 with a step of 0.001 to validate SAVI's robustness using a negative soil adjustment factor. Both MODIS and field measured LAI data were used for the validation. The results of the linear regression are shown in Figure 6.
From both LAI data, we can observe that the best results mainly appeared for X SAVI between −0.18 and 0, indicating that negative X SAVI values were very satisfactory and behaved well compared to positive values. In the MODIS LAI data (Figure 6a), the negative R 2 and slope appeared in the interval from 0.225 to 1 for X SAVI , indicating that the correlation between SAVI and LAI was extremely weak. Besides, the negative slope meant a negative relationship between LAI and SAVI. The maximin R 2 (0.2472) and minimum P (0.0003) appeared at X = −0.148. The negative R 2 and slope in the field measured data (Figure 6b) appeared for X SAVI between −0.3 and −0.2. Besides, R 2 and P were found to be highly volatile in this region. The maximin R 2 (0.6417) and minimum P (0.0863) appeared at X SAVI = −0.183. With X SAVI from −0.184 to 1, the performance of SAVI gradually worsened, with a decreasing R 2 and an increasing P. From both LAI data, we can observe that the best results mainly appeared for XSAVI between −0.18 and 0, indicating that negative XSAVI values were very satisfactory and behaved well compared to positive values. In the MODIS LAI data (Figure 6a), the negative R 2 and slope appeared in the interval from 0.225 to 1 for XSAVI, indicating that the correlation between SAVI and LAI was extremely weak. Besides, the negative slope meant a negative relationship between LAI and SAVI. The maximin R 2 (0.2472) and minimum P (0.0003) appeared at X = −0.148. The negative R 2 and slope in the field measured data (Figure 6b) appeared for XSAVI between −0.3 and −0.2. Besides, R 2 and P were found to be highly volatile in this region. The maximin R 2 (0.6417) and minimum P (0.0863) appeared at XSAVI = −0.183. With XSAVI from −0.184 to 1, the performance of SAVI gradually worsened, with a decreasing R 2 and an increasing P.  According to a earlier study [11], the optimal value of the soil adjustment factor is related to the general condition of the vegetation, with X SAVI close but smaller than one for sparse vegetation and close to but not smaller than 0 for dense vegetation. However, the lower boundary zero for dense vegetation was questioned by Ren et al. [38] based on field-measured data in arid grasslands and by Zhen et al. [8,14] based on simulated data by the discrete anisotropic radiative transfer [39][40][41][42][43][44][45][46][47][48] (DART, https://dart.omp.eu/#/, accessed on 11 March 2021) model. Remote sensing models such as DART that simulate the reflectance of land surfaces are very useful tools to investigate the domain of validity of VIs. More generally, they have a great potential to invert satellite observations in terms of LAI maps.
Similar to Ren et al. [38], we found an optimal performance of the negative soil adjustment factor. However, the vegetation conditions differed from the previous study [38] where arid grasslands gave had a NDVI of around 0.2, indicating very sparse vegetation. In our study area, the NDVI values were mainly distributed around 0.8, indicating very dense vegetation. We assumed that the better performance of negative soil adjustment factors in the two situations is because the intersections between the soil line and vegetation isolines were located in the first quadrant, even though their vegetation density varied dramatically.
This paper further confirmed that the lower boundary of the soil adjustment factor could be negative, based on field-measured and remote sensing data ( Figure 6). Also, as was already pointed out by Zhen et al. [8] based on simulated data, the determining factors of the location of intersections between the soil line and vegetation isolines are determined by the interception of the vegetation isolines and the soil lines. When the vegetation isoline intercept is larger than the soil line intercept, the intersection point is located in the second or third quadrant of the red-NIR plane. When the vegetation isoline intercept is smaller than the soil line intercept, the intersection point is located in the first quadrant. Moreover, the intercept depends also largely on the spectral properties of the leaves and soil apart from canopy structure parameters such as LAI. Our study further confirmed this conclusion using both field-measured and remote sensing data.
Besides, we noticed that NDVI had a better performance than SAVI, which is considered an improved version of NDVI. Similar undesired SAVI performances have also been reported by [38,49]. Reasons for this poor performance may contribute to the fact that NDVI assumes vegetation isolines converge to the original point, and SAVI assumes vegetation isolines converge to one point (−0.5, −0.5) located in the third quadrant. In our study (very dense vegetation areas), the vegetation isolines converged to one point in the first quadrant. As a result, the assumption based on the original point is much better than the point (−0.5, −0.5). This reason can also be seen in the TSAVI slightly better performance than SAVI because TSAVI assumes that the vegetation isolines' intersection lies between the origin and the point (−0.5, −0.5). Another possibility is that NDVI was used in the backup algorithm for inversion in MODIS LAI data when the radiative transfer inversion method fails, which might also influence the results.

Saturation Effect
Numerous studies [8,[11][12][13]15,16] have indicated that VIs show significant saturation effects with LAI growth. Several methods have been implemented to enhance the saturation resistance by adding spectral information or geometry information [19,36,50]. For example, EVI added a blue band which is more sensitive than the red band for high LAI [19]. The absorption intensity of the leaf in the blue band is less than that in the red band. When the LAI exceeded the threshold value, the canopy reflectance had almost no change in the red band due to the strong absorption intensity, and the canopy reflectance still decreased with the LAI because of the weaker absorption intensity in the blue band. Also, modified normalized difference vegetation index (MNDVI) [50] involved a combination of narrow bands in the shorter wavelengths of the red edge (700-750 nm) and longer wavelengths of the red edge (750-780 nm). Besides, NHVI [36] added an angular index representing the distribution of leaves inside the canopy to enhance the saturation resistance ability.
However, the increased information required leads to higher demands on the sensors. For example, EVI products cannot be produced from a sensor without the blue band (i.e., the advanced very high resolution radiometer, AVHRR). The narrowband requirement makes MNDVI hard to be calculated from some common sensors (i.e., Landsat 8). NHVI requires the sensor can capture multiple-angler observation. Therefore, it is a challenge for these VIs to produce long-term products across sensor systems with variable spectral response functions, swath width, and orbiting geometry. Compared with these methods, adjusting the soil adjustment factor may be effortless because it requires less information from the sensors. However, although the method using negative soil adjustment factor achieves better advantages in LAI estimation than other VIs, the accuracy is still lacking compared to these LAI estimation methods. Our method is, therefore, more suitable for use in situations where sensor conditions are limited. For example, when only an optic sensor with a fixed viewing direction is available.
In addition, we found that EVI and SAVI were negatively correlated with LAI in the MODIS LAI dataset. This negative correlation was not observed in the field measured data. It leads to inaccurate LAI estimation with errors that cannot be easily estimated. The negative correlation between VIs and the LAI is few reported in other VIs saturation studies, and we did not find a solid explanation for it. We simply assume that it is related to the insensitivity of reflectance values in the red band in the high LAI region. It will need to be further investigated in the future.

Conclusions
The saturation effects limit the use of VIs. In this study, we explored the possibility of using the negative soil adjustment factor to mitigate the saturation effects. Two data sets were used for the validation, including a long time series observation of MODIS LAI data in the Apiacás area and a global LAI dataset. The surface reflectance data derived from Sentinel-2, Landsat-8, and Landsat-5 were used for VIs calculation. An exhaustive search was employed to estimate the optimal soil adjustment factor for LAI estimation.
Results show that LAI was highly cyclical with seasonal changes. Besides, VIs calculated using the positive soil adjustment factor had serious saturation effects, and the correlation between VIs and LAI was very low. In the MODIS LAI data, SAVI and EVI even showed a negative growth with LAI. The exhaustive search suggested that a negative X SAVI (e.g., −0.148) in MODIS LAI and X SAVI (e.g., −0.183) in the field measured dataset might be better than the commonly used X SAVI (e.g., 0.5) in a high LAI area.
So far, apart from the arid grasslands and simulated data that have been mentioned, the feasibility of using negative soil adjustment factors, our study indicates the feasibility of using a negative soil adjustment factor for saturation resistance in dense vegetation areas. Our study further confirms that the lower bound of the soil adjustment factor can be less than zero and reach negative regions. However, although the accuracy of using negative soil conditioning factors is much higher than positive values, the accuracy of our estimates still falls short compared to some other methods that are more demanding on the sensor. Therefore, the method proposed in this paper is more suitable for use with common optical satellites with single fixed observation direction.