Characterizing the Error and Bias of Remotely Sensed LAI Products: An Example for Tropical and Subtropical Evergreen Forests in South China

: Leaf area is a key parameter underpinning ecosystem carbon, water and energy exchanges via photosynthesis, transpiration and absorption of radiation, from local to global scales. Satellite-based Earth Observation (EO) can provide estimates of leaf area index (LAI) with global coverage and high temporal frequency. However, the error and bias contained within these EO products and their variation in time and across spatial resolutions remain poorly understood. Here, we used nearly 8000 in situ measurements of LAI from six forest environments in southern China to evaluate the magnitude, uncertainty, and dynamics of three widely used EO LAI products. The ﬁner spatial resolution GEOV3 PROBA-V 300 m LAI product best estimates the observed LAI from a multi-site dataset ( R 2 = 0.45, bias = − 0.54 m 2 m − 2 , RMSE = 1.21 m 2 m − 2 ) and importantly captures canopy dynamics well, including the amplitude and phase. The GEOV2 PROBA-V 1 km LAI product performed the next best ( R 2 = 0.36, bias = − 2.04 m 2 m − 2 , RMSE = 2.32 m 2 m − 2 ) followed by MODIS 500 m LAI ( R 2 = 0.20, bias = − 1.47 m 2 m − 2 , RMSE = 2.29 m 2 m − 2 ). The MODIS 500 m product did not capture the temporal dynamics observed in situ across southern China. The uncertainties estimated by each of the EO products are substantially smaller (3–5 times) than the observed bias for EO products against in situ measurements. Thus, reported product uncertainties are substantially underestimated and do not fully account for their total uncertainty. Overall, our analysis indicates that both the retrieval algorithm and spatial resolution play an important role in accurately estimating LAI for the dense canopy forests in Southern China. When constraining models of the carbon cycle and other ecosystem processes are run, studies should assume that current EO product LAI uncertainty estimates underestimate their true uncertainty value.


Introduction
The Leaf Area Index (LAI), the total leaf area per unit ground area, is a key biophysical variable playing an important role in global carbon, water, and energy cycles [1,2].As such, it acts as an important parameter for several applications, such as land surface models [3], ecological models [4], and yield prediction models [5].The amount of leaf area has a first-order control on photosynthesis, transpiration, and absorption of radiation, varying in both space and time.LAI seasonal dynamics provide information about phenological processes of canopy development, senescence, and plant traits [6].Therefore, retrieving LAI over large areas and having a good knowledge of their yearly variations, errors, and bias is extremely important.Such information is central to accurately estimating primary productivity, understanding land surface-atmosphere exchanges, and detecting the response of terrestrial vegetation to climate change [7].It is also beneficial for a large remote sensing community because it provides insights for the interpretation and correct usage of the LAI maps.Satellite-based Earth Observation (EO) offers the opportunity to retrieve information on LAI which is global in coverage and at increasing temporal and spatial resolution [8][9][10][11][12].Robust estimates of uncertainties associated wth EO LAI estimates are needed to ensure their appropriate integration within further analyses such as data assimilation [6,13].However, such robust evaluations are rarely possible due to the scarcity of in situ data at appropriate temporal and spatial resolutions [14,15].
LAI product uncertainty information represents the performance of the products' algorithm and reflects the uncertainties in the input data, model imperfections, and the inversion process [16][17][18], which could be called the theoretical uncertainties [19].Evaluation studies typically assess the theoretical uncertainty of LAI based on the standard from Global Climate Observing System (GCOS) [20], which sets a target for absolute LAI uncertainty range to be within 0.5 m 2 m −2 and a relative uncertainty of less than 20%.However, few studies validate the uncertainty estimates of EO LAI using in situ data [14] at appropriate scales.A key challenge for the validation process is the mismatch in the scale of satellite products (typically > 300 m) compared to that of the in situ data (<10 m).Heterogeneity of LAI at scales finer than the satellite resolution [21] means that simple comparison between measurements at different scales can be problematic [22].
The validation of LAI at temporal scales includes a more important component, the LAI temporal dynamics, besides the traditional absolute value checking [23].The temporal dynamics of LAI time series contain useful ecological information (e.g., phenological infomation), which help study the relationship between plant phenology and climate [24] or assist in land cover classification [25].However, the temporal dimension has been neglected for most of the LAI validation work [23].To validate the seasonal variations of LAI products, higher field sampling rates are need, especially for the Evergreen Broadleaf Forests (EBF) across several regions (Africa, Eurasia, South America, Australia and Asia) [14].
Despite the important role played by tropical forests in regulating the global carbon [26,27] and energy balance [28], tropical phenology is highly uncertain.Large-scale monitoring of tropical forest LAI with satellites is hindered by several challenges, including signal saturation [23], poor observation conditions (cloud-aerosol contamination, etc.) [29], coarse spatial and temporal resolutions [30], and scarcity of both ground data and detailed validation tests [31][32][33].Uncertainties are therefore particularly high in these regions [34].Some progress has been made in the Amazon [35][36][37][38][39] and African forests [40][41][42]; however, the Southern China region has typically been neglected in remote sensing phenology studies [43][44][45][46].Consequently, the phenological character of forests in this region remains uncertain, with additional complexity driven by fragmentation [47], high species diversity, and complex topography [48].The lack of detailed validation studies in Southern China means that there is little information regarding the suitability of global EO LAI products for the forests in this region, which play a major role in the carbon sink across China, and span the climate gradient from the subtropics to tropics [49,50].
In this study, we take advantage of a network of 6 sites with ~8000 ground-based measurements of LAI for forests located across the tropical and subtropical Southern China region to fully evaluate three satellite-based remote sensed LAI products with global coverage (Moderate Resolution Imaging Spectroradiometer (MODIS) 500 m and PROBA-V GEOV2 1 km and V3 300 m).Specifically, we address the following questions: (i) How well do EO-based estimates of LAI capture temporal dynamics observed in ground measurements?(ii) How robust are EO LAI error estimates?
To answer these questions, we evaluate errors across multiple LAI retrieval algorithms and spatial resolutions and their associated temporal dynamics.We use additional fine resolution satellite data to evaluate the heterogeneity of canopy cover at the scales intermediate between in situ data and the satellite products at each location.Finally, we discuss the importance of robust error characterization for LAI products in the context of constraining models of the terrestrial carbon cycle and other ecosystem processes.We aim to characterize and understand errors and bias in the EO LAI products but do not make retrievals over the entire region.Whilst our scientific questions are focused on the Southern China region, the approach used can be applied elsewhere if in situ data is available.Assessing the uncertainties in LAI products through comparison with in situ measurements, i.e., direct validation is critical for their proper use in land surface models [51,52].A better understanding of the uncertainties embedded in current LAI products will improve the assimilation of the LAI into land surface modeling studies [53,54] by providing a more robust error weighting [55].

Study Area
Site-level LAI observational data are collected from CERN (Chinese Ecosystem Research Network) [56].The six sites included are Ban Na Forest (BNF), Ai Lao Forest (ALF), Gong Ga Forest (GGF), Ding Hu Forest (DHF), He Shan Forest (HSF), and Shen Nong Forest (SNF).These sites are spread across southern China between 21.9 • and 31.3 • N, covering major plant functional types from northern subtropical zones to the northern tropical zones in this region (Figure 1).The regional climate is characterized by a wet, warm summer and a dry, mild winter [57].The site-specific Koeppen climate classifications [58] are Cwa (Temperate zone warm summer; dry winter) for ALF, BNF, HSF, and SNF.Cfa (Temperate zone hot summer; no dry season) for DHF and Dwb (Cold zone dry winter; cold summer) for GGF.The mean annual temperature ranges between 4.2 and 21.8 • C. Annual precipitation ranges between 1506 mm and 2175 mm.The sites vary in elevation (above sea level) between 70 and 3160 m (Table 1).Forest types include the tropical seasonal rainforests, mixed coniferous forest, evergreen and deciduous mixed broad-leaved forest.The sites used support publically forests which are natural, except for forests at HSF and DHF, which support planted forests over 40 years old.All forests have at least three vertical layers: woody, shrub, and herbs.The tree density (trees per hectare) varies widely across sites.The highest tree density reaches over 7000 and the lowest is around 600 (Table S1).The range of the mean diameter breath heights is 5-20 cm (Table S1).The dominant species composition differs between sites but mainly consists of the evergreen species (Table S1).An image of forests measured at each site is included in Figure S1. the satellite products at each location.Finally, we discuss the importance of robust error characterization for LAI products in the context of constraining models of the terrestrial carbon cycle and other ecosystem processes.We aim to characterize and understand errors and bias in the EO LAI products but do not make retrievals over the entire region.Whilst our scientific questions are focused on the Southern China region, the approach used can be applied elsewhere if in situ data is available.
Assessing the uncertainties in LAI products through comparison with in situ measurements, i.e., direct validation is critical for their proper use in land surface models [51,52].A better understanding of the uncertainties embedded in current LAI products will improve the assimilation of the LAI into land surface modeling studies [53,54] by providing a more robust error weighting [55].

Study Area
Site-level LAI observational data are collected from CERN (Chinese Ecosystem Research Network) [56].The six sites included are Ban Na Forest (BNF), Ai Lao Forest (ALF), Gong Ga Forest (GGF), Ding Hu Forest (DHF), He Shan Forest (HSF), and Shen Nong Forest (SNF).These sites are spread across southern China between 21.9° and 31.3°N,covering major plant functional types from northern subtropical zones to the northern tropical zones in this region (Figure 1).The regional climate is characterized by a wet, warm summer and a dry, mild winter [57].The site-specific Koeppen climate classifications [58] are Cwa (Temperate zone warm summer; dry winter) for ALF, BNF, HSF, and SNF.Cfa (Temperate zone hot summer; no dry season) for DHF and Dwb (Cold zone dry winter; cold summer) for GGF.The mean annual temperature ranges between 4.2 and 21.8 °C .Annual precipitation ranges between 1506 mm and 2175 mm.The sites vary in elevation (above sea level) between 70 and 3160 m (Table 1).Forest types include the tropical seasonal rainforests, mixed coniferous forest, evergreen and deciduous mixed broad-leaved forest.The sites used support publically forests which are natural, except for forests at HSF and DHF, which support planted forests over 40 years old.All forests have at least three vertical layers: woody, shrub, and herbs.The tree density (trees per hectare) varies widely across sites.The highest tree density reaches over 7000 and the lowest is around 600 (Table S1).The range of the mean diameter breath heights is 5-20 cm (Table S1).The dominant species composition differs between sites but mainly consists of the evergreen species (Table S1).An image of forests measured at each site is included in Figure S1.

Field LAI Data Measurements
In each of the six sites, there are one to four subplots at which LAI is estimated, and each subplot area ranged from 400-10,000 m 2 (Figure S1).Field measurements of LAI were made each month using an LAI-2000 [59].Field LAI measurement processing is consistent for each site and is based on the unified data collection and quality control protocols specified for CERN [60].The sampling positions are distributed evenly and diagonally within the plot.The horizontal distance between each sampling points or plot boundary are to exceed 15 m and 10 m separately to avoid overlap sampling and reduce marginal effects.LAI-2000 were used to scan the canopy twice a day (8:00 am and 16:00 pm) to get the three-layer LAI (woody, shrub, and herbs) at each sampling point at the measurement day.The position for each sampling measurements and scan height is fixed.The median number of monthly measurements made at each site ranged from 6 to 80 (Table S2).BNF is the most intensively surveyed site, with between 203 and 360 measurements per year; therefore, BNF is particularly suited to evaluating the temporal dynamics of EO LAI products.For the other sites, only a subset of years and growth season months are measured (Table S2).A total of 795, 9 ground measurements from 2005 to 2017 were acquired across all the sites (Table S2).

Earth Observation LAI Estimates
In this study we evaluate three global EO LAI products.MODIS 500m product (MOD15A2H) provides estimates of LAI at 500 m spatial resolution with an 8-day interval from 2000-present [10].GEOV2 1 km [11] and GEOV3 300 m LAI [12] were derived from the SPOT/VEGETATION sensor data at a 10-day interval, at spatial resolutions of 1/112 • and 1/336 • , respectively (approximately 1-km and 300 m at the equator), in the Plate Carrée projection.For the MODIS 500 m LAI product, we consider only the products derived from the main algorithm, which is based on the use of Look Up Tables (LUTs) built for six different plant functional types, with simulations from a three-dimensional radiative transfer model [61].For the GEOV2 1 km and GEOV3 300 m LAI products, LAI is estimated using Neural Networks(NNTs) applied on Top-of-Canopy (TOC) input reflectance in the red, near-infrared, and shortwave infrared bands, at 1km resolution and 300m, respectively [11,12].The period of MODIS 500 m LAI products and GEOV2 1 km LAI products is from 2005 to 2017 and 2014 to 2017 for the GEOV3 300 m LAI product (released from 2014).
The definition of uncertainty information varies between LAI products.For MODIS LAI, the standard deviation is used to measure the uncertainty of pixel LAI values, which is calculated over all acceptable solutions of a look-up table (LUT) retrieval method [17].For GEOV2 1 km and GEOV3 300 m LAI products, the uncertainties are computed as the RMSE between the final decadal value and the daily NNTs estimates in the compositing period [12,62].In this study, they are both treated as the LAI products' uncertainty.The quality assessment information is used to clean the LAI and inform uncertainty information at the pixel level for all of these LAI products.For MODIS 500 m products, data effected by cloud shadow, internal cloud masks, or aerosol are removed and only main retrievals were used in the analysis.For Copernicus LAI data (GEOV2 1 km and GEOV3 300 m), pixels that were filled or interpolated, or out of LAI range, were removed.

Mapping Heterogeneity of Canopy Properties
To evaluate the upscaling of field LAI measurements to satellite products, Landsat Level-2 Surface Reflectance products including Enhanced Vegetation Index (EVI) and Normalized Difference Vegetation Index (NDVI) derived from Landsat 5 (TM), Landsat 7 (ETM+) [63] and Landsat 8 (OLI) [64] scenes are also analyzed.EVI is calculated from red (R), near-infrared (NIR) and blue band (B) based on Equation ( 1) and incorporates parameters to adjust for canopy background, atmospheric resistance.
NDVI is calculated as a ratio between the R and NIR values in a traditional fashion (Equation ( 2)).
Landsat 5 (TM), Landsat 7 (ETM+), and Landsat 8 (OLI) Level-2 surface reflectance products were generated, using the Land Surface Reflectance Code (LaSRC) and the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithms [65,66].The LEDAPS and LaSRC surface reflectance algorithms correct for the temporally, spatially, and spectrally varying scattering and absorbing effects of atmospheric gases and aerosols, which is essential to derive the Earth's reflectance surface values.Landsat Spectral Indices are generated at 30 m spatial resolution on a Universal Transverse Mercator (UTM) mapping grid.The Spectral Indices could be obtained from the USGS Earth Resources Observation and Science (EROS) on-demand processing system.
The GTOPO 30 1 km elevation map [67] in the China region with Landsat scene footprints at each site is presented in Figure S2.In each of the six sites, there are one to four sampling subplots.The location of each site with the sampling subplots is shown on the Landsat scenes with path and row information (Figure S3) and on the 90 m SARTM elevation map [68], respectively (Figure S4).There are in total 778 Landsat 5-7 and Landsat 8 scenes with good quality (collection category is Tier 1) and cloud cover smaller than 20% of the total scene area from the year of 2005 to 2017 for all of these forests.Because there is no quality assessment information for the Landsat spectral indices, we use the pixel quality assurance band (pixel_qa) instead to obtain quality information.Pixels with cloud shadow, aerosols, etc., are filtered out, and we keep pixels with the best quality.
Landsat EVI and NDVI were extracted at the field sampling plots spatial level (20~100 m).Based on the Landsat scene quality information, we cleaned and kept just good quality Landsat pixel data.Then, daily EVI and NDVI were retrieved from the area-weighted averaged Landsat pixels inside the field sampling plot.We averaged the daily EVI and NDVI values by month.Finally, the monthly EVI and NDVI at the field sampling plots level were regressed with corresponding monthly in situ LAI values to calculate the statistics (R 2 , etc.).
We investigated whether forest cover fraction in different LAI products' spatial resolution levels influenced product bias and error.We use the 30 m resolution GlobeLand30 land cover product for 2010 [69] to extract the forest cover fraction inside the relevant pixels for each of the LAI products at each site.We classify the land cover for each 30 m pixel on different spatial scales (LAI products' scale: 300 m, 500 m, and 1 km) at each site and plot.We only kept the 30 m pixels located entirely inside the LAI pixels.The total number of 30 m pixels and pixels classified as forests were both counted and summed to determine the total area and the forest areas at different spatial scales.Then we can obtain the forests cover proportion at different spatial scales for each plot and site.Finally, the forests' cover proportion was compared with the LAI products' bias at each site.

Statistical Analysis
Daily satellite LAI and uncertainty information were retrieved from the three LAI products at the location of each subplot of the six sites (Figure S1).If the field subplots overlapped with multiple pixels for a given LAI product, then the satellite-based value was estimated using the weighted sum of the contributing pixels, weighted by the area of overlap (Equation ( 3)).Data were cleaned based on the products' quality assessment information, before aggregation into monthly mean values, and comparison with the mean monthly ground-based LAI measurements at each site and plot.We calculate R 2 , bias, and RMSE between the field and satellite LAI time series to evaluate the accuracy and error of the three LAI products.
where LAI R is the extracted remote sensing LAI value; n represents the total number of the contributing pixels which overlapped with the field plots; A is the total area for the plots, Ai is the overlapped area between pixel and plots; Pi is the remote sensing pixel value.LAI products' uncertainty information was evaluated using the ratio of LAI products' uncertainty (LAI U ) with the bias against the field measurements, r, as follows (Equation (4)), where the bias (Equation ( 5)) is defined as: The ratio, r, can be used to evaluate whether or not the product uncertainty is robust.When r ≥ 1, the bias in the LAI product is within the reported uncertainty range, indicating that the uncertainty estimate for LAI is robust; conversely, if r < 1, the reported uncertainty fails to capture the observed bias and therefore the uncertainty estimate is not robust.We calculate this ratio r for each month at different plots and sites.

Calculation of Amplitude, Phases and Periods for LAI Time Series
We use metaCycle [70] to calculate the periodic characteristics (periods, phases, and amplitude) for the monthly field, MODIS, GEOV2 and V3 LAI time series.This method can calculate the periodic characteristics for time series with missing values [70], which frequently occurred for the satellite LAI times series.MetaCycle incorporates three methods: ARSER (Autoregressive spectral analysis) [71], JTK_CYCLE(Jonckheere-Terpstra-Kendall algorithm) [72], and Lomb-Scargle [73], for periodic signal detection, and it could output integrated analysis.The integrated period from MetaCycle is an arithmetic mean value of multiple periods, while phase integration based on the mean of circular quantities calculated from the above three mentioned methods.
MetaCycle recalculates the amplitude with the following model (Equation ( 6)): where Yi is the observed value at time t i ; B represents the baseline value (mean value) of the LAI time series; TRE is the trend level of the time-series profile; A is the amplitude of the waveform.PER and PHA are the integrated period and phase, respectively.PER represented the periods of the LAI time series; here, it is 12 months.PHA represented the phases of the LAI time series at the period of 12 months."i" is the time points of the time series."n" represents the total number of the time points of the time series.In this model, only B, TRE, and A are unknown parameters and can be calculated using ordinary least squares.Fisher's method is implemented in MetaCycle for integrating multiple p-values.
Relative amplitude (rAMP), which is the A/B, can be used to compare the amplitudes between time series with different baselines levels.
In this study, each LAI time series is decomposed into the baseline value (mean value), the trend, and the wavelength with specific amplitude, period, and phase.The integrated phases can be used to compare whether the LAI time series is lead, lag, or synchronous.The amplitude can be used to quantify the magnitude of seasonal fluctuations for LAI time series and the relative amplitude can be used to compare the LAI seasonal fluctuations from different sites and data source.We apply this method to the monthly averaged field LAI time series at BNF, and the monthly averaged retrieved satellite LAI at each site because only BNF has frequent enough field measurements over multiple months from the year 2005 to 2017.

Software
Remote-sensing data processing was carried out by Rv3.4.0 [74] using the package "raster" [75].Amplitude, phases, and periods of LAI time series were calculated using the functions "meta2d" in the package 'MetaCycle' [70] in R. Figures and maps were produced in R and ArcGIS 10.4.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 20 used to compare the LAI seasonal fluctuations from different sites and data source.We apply this method to the monthly averaged field LAI time series at BNF, and the monthly averaged retrieved satellite LAI at each site because only BNF has frequent enough field measurements over multiple months from the year 2005 to 2017.

Software
Remote-sensing data processing was carried out by Rv3.4.0 [74] using the package "raster" [75].Amplitude, phases, and periods of LAI time series were calculated using the functions "meta2d" in the package 'MetaCycle' [70] in R. Figures and maps were produced in R and ArcGIS 10.4.

Accuracy of the LAI Products Magnitude
Satellite derived mean LAI estimates were typically lower than the field measurements, which varied between 3.2 and 5.6 m 2 m −2 (Table 2).GEOV3 300 m estimates were 10% lower than field measurements (mean = 4.50 m 2 m −2 ), MODIS 500m were 30% lower (mean = 3.6 m 2 m −2 and GEOV2 1 km were 40% lower (mean = 3.07 m 2 m −2 ) (Figure 2; Table 2).LAI mean and relative uncertainty increased with product resolution: GEOV3 300m (0.22, 7%) < MODIS 500 m (0.37, 10%) < GEOV2 1 km (0.56, 19%).m −2 ) and their statistical metrics calculated using months which have both field and EO LAI values from 2005 to 2017."r" is median of the ratio between product uncertainty and bias for all months and subplots.R 2 was calculated from the linear regression between satellite products and field monthly measurements."N" denotes the number of monthly values for all subplots at each site.All the statistics are calculated on a monthly scale.Field mean LAI was best captured by GEOV3 300 m (R 2 = 0.45; Figure 3), followed by GEOV2 1 km (R 2 = 0.36), and MODIS 500 m (R 2 = 0.20) for all of these six forests in this region.This performance varied between forests and all of these products showed the highest R 2 (Table 2) for the forests at SNF, which showed more LAI seasonal variation (Figure 2).GEOV3 300 m LAI also had the lowest RMSE and bias (RMSE = 1.21, bias = −0.41)compared with GEOV2 1 km LAI (RMSE = 2.32, bias = −2.04)and MODIS 500m LAI (RMSE = 2.29, bias = −1.47)(Figure 3, Table 2).

Robustness of the LAI Products Uncertainty
The magnitude of the reported EO LAI uncertainty is typically smaller than the bias to in situ estimates (Figure 4), here shown as the ratio of uncertainty and reported bias (r), and this result is consistent across forests (Figure S5, Table 2).GEOV3 300 m has the smallest proportion of EO uncertainty estimates below the observed bias (86%) while MODIS 500m LAI and GEOV2 1 km have 87% and 91% of uncertainty estimates below the reported bias.The median ratio indicated similar average skill between both GEOV3 300 m (0.26) and GEOV2 1km LAI products (0.28) and MODIS 500 m LAI products (0.21) (Figure 4, Table 2).

Robustness of the LAI Products Uncertainty
The magnitude of the reported EO LAI uncertainty is typically smaller than the bias to in situ estimates (Figure 4), here shown as the ratio of uncertainty and reported bias (r), and this result is consistent across forests (Figure S5, Table 2).GEOV3 300 m has the smallest proportion of EO uncertainty estimates below the observed bias (86%) while MODIS 500m LAI and GEOV2 1 km have 87% and 91% of uncertainty estimates below the reported bias.The median ratio indicated similar average skill between both GEOV3 300 m (0.26) and GEOV2 1km LAI products (0.28) and MODIS 500 m LAI products (0.21) (Figure 4, Table 2).

Validation of the LAI Products Temporal Dynamics
The metaCycle analysis showed that a statistically significant period of 12 months (p < 0.05) is present in the monthly field LAI time series at BNF and in all monthly satellite LAI time series at all sites, except for MODIS 500 m LAI, which showed no significant 12-month periodic characteristics at BNF.Overall, LAI retrieved from GEOV3 300 m showed very close values of the base, phase, amplitude and relative amplitude with the field LAI at BNF (Figure 5).LAI retrieved from GEOV2 1km showed a lower base and amplitude, similar phase, but larger relative amplitude compared with GEOV3 300 m LAI for most sites and the field at BNF (Figure 5).LAI retrieved from MODIS 500 m product showed the lower base values, but had a high variation of phases with one underestimated phase for the forest at ALF compared with the other two LAI products (Figure 5b).The relative amplitude of LAI retrieved from MODIS 500 m is lower for the forest at GGF and SNF, but higher at HSF compared with the other two LAI products (Figure 5c).

Validation of the LAI Products Temporal Dynamics
The metaCycle analysis showed that a statistically significant period of 12 months (p < 0.05) is present in the monthly field LAI time series at BNF and in all monthly satellite LAI time series at all sites, except for MODIS 500 m LAI, which showed no significant 12-month periodic characteristics at BNF.Overall, LAI retrieved from GEOV3 300 m showed very close values of the base, phase, amplitude and relative amplitude with the field LAI at BNF (Figure 5).LAI retrieved from GEOV2 1km showed a lower base and amplitude, similar phase, but larger relative amplitude compared with GEOV3 300 m LAI for most sites and the field at BNF (Figure 5).LAI retrieved from MODIS 500 m product showed the lower base values, but had a high variation of phases with one underestimated phase for the forest at ALF compared with the other two LAI products (Figure 5b).The relative amplitude of LAI retrieved from MODIS 500 m is lower for the forest at GGF and SNF, but higher at HSF compared with the other two LAI products (Figure 5c).

Evaluation of Landcover Heterogeneity
The Landsat data provide a more spatially consistent comparison of surface reflectance against field data of LAI.The results (Figure S6) show that there are significant correlations between NDVI/EVI and LAI across sites, but the effect sizes are small (R 2 = 0.1).Thus, the relationships between remotely-sensed NDVI/EVI at the fine resolution and field-measured LAI are less significant than those between LAI products at the coarse resolution and field-measured LAI.The forest cover proportion at the pixel level varied between 38%-100% (Table S3).The landcover data show that LAI retrievals were degraded by the heterogeneity of landcover within product pixels.There was a clear pattern of LAI bias increasing with reductions in the proportions of forest cover at the pixel scale (Figure S7).

Discussion
The performance of three EO LAI products was evaluated at six forests across southern China.GEOV3 300 m LAI showed the best fit (R 2 = 0.45) and the smallest bias (bias = −0.54)(Figure 3) and the most similar LAI time series seasonal dynamic variation (phase and amplitude) compared with in situ LAI measurements (Figure 5).For the GEOV2 1 km LAI product, its performance had larger bias (Figure 3) and seasonality but similar phases of the LAI time series with the in situ observations compared with GEOV3 300 m LAI products (Figures 3 and 5).
The different performance of the two Copernicus products is likely due to the mismatch of the scaling between the EO pixel size and the field site.LAI is strongly non-linearly related to reflectance, making its estimation from remote sensing observations scale-dependent [76,77].In contrast, the core operational algorithm (neural network techniques), data filtering and smoothing processes are similar for these two products.There are differences in the method used for temporal compositing, where temporal smoothing and gap filling using a climatology are used for GEOV2 1km and A detailed description of the calculation can be found in the Methods section of this paper.

Evaluation of Landcover Heterogeneity
The Landsat data provide a more spatially consistent comparison of surface reflectance against field data of LAI.The results (Figure S6) show that there are significant correlations between NDVI/EVI and LAI across sites, but the effect sizes are small (R 2 = 0.1).Thus, the relationships between remotely-sensed NDVI/EVI at the fine resolution and field-measured LAI are less significant than those between LAI products at the coarse resolution and field-measured LAI.The forest cover proportion at the pixel level varied between 38%-100% (Table S3).The landcover data show that LAI retrievals were degraded by the heterogeneity of landcover within product pixels.There was a clear pattern of LAI bias increasing with reductions in the proportions of forest cover at the pixel scale (Figure S7).

Discussion
The performance of three EO LAI products was evaluated at six forests across southern China.GEOV3 300 m LAI showed the best fit (R 2 = 0.45) and the smallest bias (bias = −0.54)(Figure 3) and the most similar LAI time series seasonal dynamic variation (phase and amplitude) compared with in situ LAI measurements (Figure 5).For the GEOV2 1 km LAI product, its performance had larger bias (Figure 3) and seasonality but similar phases of the LAI time series with the in situ observations compared with GEOV3 300 m LAI products (Figures 3 and 5).
The different performance of the two Copernicus products is likely due to the mismatch of the scaling between the EO pixel size and the field site.LAI is strongly non-linearly related to reflectance, making its estimation from remote sensing observations scale-dependent [76,77].In contrast, the core operational algorithm (neural network techniques), data filtering and smoothing processes are similar for these two products.There are differences in the method used for temporal compositing, where temporal smoothing and gap filling using a climatology are used for GEOV2 1 km and interpolation applied in GEOV3 300 m) [78].Differences in the applied gap-filling approach between these two products do not impact our conclusion that resolution is the primary driver of performance improvement at 300 m relative to 1 km, as all gap-filled and interpolated retrievals were removed in our study.The in situ measurements were usually conducted in protected mature forest areas with high plant densities, complex canopy stratums, and rich species diversities (Figure S1, Table S1).However, at coarser spatial resolutions, pixels integrate heterogeneity over a greater diversity of landscapes, habitats, and species, distributed across a variety of stages of growth and succession [79].Thus, fine resolution GEOV3 300 m LAI product showed lower bias and similar seasonal dynamic variation to the in-situ measurements compared to the higher bias for the GEOV2 1 km LAI product.
For Collection 6 500 m MODIS LAI, our results indicate that this product does not capture the dynamics observed in situ, irrespective of the accuracy of estimated in situ LAI (Figure 3).Furthermore, it does not adequately capture the LAI seasonal dynamics (Figure 5).This result is consistent with two validation studies which both showed MODIS LAI had the poorest performance for the evergreen forests in the south of China [80,81].The MODIS LAI values for tropical evergreen forests are severely impacted by atmospheric conditions, especially clouds during the growing season (around 42% data are influenced by the cloud in this study, Figure S8), which lead to strong noise in the input reflectance data and affect the retrieval [78].Additionally, the reflectance saturation usually happened in dense canopies and the main algorithm is sensitive to uncertainties in atmospheric correction, particularly when red and NIR BRFs are saturated [82,83].This means the reflectance does not contain sufficient information to estimate the LAI value [84] and leads to the instability of the LAI retrievals [14].All of these may result in the poor performance of representation of evergreen forests in southern China for the MODIS 500 m LAI products.Furthermore, the product does not adequately capture LAI periodic characteristics in comparison to the field data and showed temporal inconsistency with GEOV2 and V3 products (Figure 5).Jiang (2017) also found the large temporal inconsistency between existing global LAI products at a longer time scale [85].Cammalleri (2019) found GEOV2 fAPAR showed a systematic overestimation of the fAPAR anomalies compared with the MODIS fAPAR and proposed a two-step harmonization procedure to remove this discrepancy [86].However, the homogenization may alter the magnitude of the original fAPAR time series in an undesirable way.These results highlight the need to validate the temporal consistency between different satellite products and explore more solutions to deal with such inconsistencies.In addition, geolocation uncertainty due to the spatial mismatch could also have influenced validation reliability, As our study sites are not homogeneous and sampling area is consistently smaller than the pixel of remote sensing data, mean or median remotely-sensed LAI values of surrounding 3 × 3 array of pixels [15,87,88] cannot be used for validation.To solve this problem, future field measurements should be conducted at a larger spatial scale and across more homogenous habitats.
Overall, the absolute and relative uncertainty of LAI products tends to be smaller for the fine resolution LAI products (Table 2).Collection 300 m GEOV3 has the smallest absolute uncertainty at around 0.22.Collection1 km GEOV2 has the largest absolute uncertainty at around 0.56.Uncertainty magnitude for Collection 6 500 m MODIS LAI is 0.37.This performance is consistent with the LAI relative uncertainty (Table 2).If compared with the absolute uncertainty requirements (±0.5) and relative uncertainty requirements (20%) set by the GCOS [12], all three products (Collection 300 m GEOV3, 0.22, 5%; Collection 6 500 m MODIS LAI, 0.37, 10%; Collection 1 km GEOV2 0.56, 18%) appear satisfactory.This is consistent with two global studies comparing different EO LAI products (MODIS, CYCLOPES, and GLOBCARBON) at two coarse spatial resolutions (5 km and 1 km).More pixels can meet the absolute uncertainty and relative uncertainty requirements at the spatial resolution of 1 km compared with the 5 km [19,34].All of these results indicated that the LAI product uncertainty is scale-dependent.The algorithms will produce smaller uncertainty estimated for the LAI values at finer spatial scales.This could make the uncertainty estimates of the fine-resolution products become very conservative.
Changes in LAI error with spatial resolution could be due to vegetation heterogeneity, especially for the forests and scale-dependent reflectance values.However, it is still unclear how product uncertainty changes over different spatial scales.For MODIS, the product uncertainty is the standard deviation over all acceptable solutions of a look-up table (LUT) retrieval method [17]; for Collection 300 m GEOV3 and 1 km GEOV2 LAI products, the quantitative uncertainties (LAI_ERR) are computed as the RMSE between the final dekadal value and the daily NNTs estimates in the compositing period [12,62].Standardizing different LAI products to the same spatial scale and averaging the LAI and uncertainty for different pixels during the upscaling and resampling process [19,89] could all generate errors [90].When dealing with different spatial scales, the land surface properties (e.g., land cover proportion, soil properties) could change significantly.Thus, LAI retrieval algorithms are based on inherent empirical assumptions on the distribution of their parameters, which can depart significantly from the actual canopy and soil characteristics [14].
The ratio (r, in Equation ( 4)) of LAI product uncertainty to the bias of in situ measurements is dramatically lower than the reference value for all of these LAI products (Figure 4, Table 2).This result indicated that none of the three LAI uncertainty products tested here were able to provide a robust estimate for the in-situ measurements.The most promising product is the Collection 300 m GEOV3 LAI; however, while its accuracy is the best, its uncertainty is still too conservative to capture the offset from co-located field LAI measurements.This weakness reflects the complexity of this biome type and highlights that the algorithm used to generate the uncertainty information for these LAI products needs to be improved.For areas of high uncertainties, data producers may need to refine the algorithms and verify the information using additional data sources, particularly in situ data [34].However, with respect to field data collection, more accurate measurements are also needed to test the results, such as destructive methods or using the litterfall traps.Direct measurement methods obtaining 'true' LAI values could also act as a reference to correct the bias and errors in indirect measurement techniques [31].
The accuracy of LAI data retrieved from remote sensing products is very important for a range of earth system models, for which LAI is a key internal variable.LAI uncertainty is particularly critical for model-data fusion in terrestrial carbon cycle studies [91].The uncertainty produced for the existing remote sensing LAI products in the analysis would further lead to biases in vegetation productivity estimates [92][93][94].For instance, the saturation of satellite-derived VIs generally results in an underestimation of GPP or NPP in areas with dense vegetation.Using finer resolution remote sensing data in carbon modeling studies could be an effective way to reduce the uncertainty in the estimates of C fluxes and/or stocks, but more robust errors are critical [95].
Uncertainty estimates of the observations are as important as the observations themselves because together they co-determine the outcome of the assimilation systems [96].Bayesian approaches are commonly used data assimilation systems and the influence of observations are weighted by observational uncertainties.Uncertainties in data are critical as they often determine the outcome of analyses and forecasts [97].The assimilation process requires clear error quantification for LAI to resolve model results and limit biases [98].Besides, the presence of bias in a data stream will limit the utility of using multiple observation types in an assimilation framework.Therefore, it is imperative to characterize the error in the observations and understand better the error associated with the direct measurements of LAI to landscape pixel [99].Our results indicated that none of these LAI uncertainty products are robust and therefore users should consider inflating the existing LAI product uncertainty when using these data sources in data assimilation frameworks [98,100,101] or in other studies for which the uncertainty of the product plays an important role.
One way of solving this spatial mismatch problem is to upscale the field measurements based on its relationship with the plant index retrieved from high-resolution images (e.g., 10-30 m resolution) [102], or using airborne Light Detection and Ranging (LiDAR) [103,104].However, we failed to find a good relationship between the Landsat based EVI or NDVI and the field measurements (Figure S5).LAI-VI relationships are limited by a number of factors, including vegetation type, sun-surface-sensor geometry, leaf chlorophyll content, background reflectance, and atmospheric quality [31].Geolocation uncertainties due to spatial mismatch and insufficient overlap data between in situ measurements and good quality Landsat scenes partly lead to this result.Saturation effects occur for both field measurement and spectral plant indices for very dense canopies, which also increases the challenge of upscaling field LAI to generate high-resolution LAI reference maps.In addition, the temporal mismatch between field and Landsat data also contributes to the poor relationship between field and Landsat data.This result may be also associated with the complex forest environments (such as terrain, tree densities, canopy structures, and species diversities) in southern China (Figure S4, Table S1).Thus, their LAI cannot be accurately modeled based on simple spectral indices, even with relatively high resolution (30 m) satellite data.More complex upscaling approaches or those utilizing very high-resolution sensors (e.g., LIDAR) to bridge the scaling gap may yield better results [102,105,106].In particular, there is an exciting potential for unmanned aerial vehicles to provide critical upscaling information from field to satellite [21].
When using EO LAI products in the absence of reliable validation data, we recommend that the provided uncertainties are treated cautiously.Remote-sensing LAI products require improved algorithms with particular attention given to generating robust uncertainty estimates.New cloud and aerosol detection techniques based on time series and spatial analysis may help to improve the uncertainty products and LAI estimated for evergreen forests [107].In addition, the new remote-sensing systems such as Global Ecosystem Dynamics Investigation (GEDI) [108], National Aeronautics and Space Administration's (NASA's) Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) and Synthetic aperture radar (SAR) provide promising ways to solve cloud coverage issues [109].Given the importance of data uncertainties (e.g., integration within terrestrial carbon modeling frameworks [92][93][94]), there is an urgent need for robust uncertainty characterization, based on links to in situ observations.The findings in this study are constrained to southern China.Future studies could explore the LAI and error products for other regions, or focus on the seasons where products show the largest uncertainties [19].

Conclusions
In this study, we validated three global LAI products and their uncertainty products using ~8000 in situ measurements of LAI from subtropical and tropical forests in the southern China region.The finest resolution product, Collection 300 m GEOV3, performed best, with the lowest RMSE and the lowest bias.It also best captured the temporal dynamics observed in the in-situ dataset, including the magnitude, amplitude, and phases of the LAI time series.Collection 6 500 m MODIS LAI showed the poorest performance and importantly it did not capture the temporal dynamics observed in situ for the forests in this specific region.Critically, for all three LAI products, the accompanying uncertainties were all far smaller than the bias compared to the in situ measurements, indicating that each product uncertainty estimate is not robust.Given the importance of LAI uncertainty in the context of constraining models of carbon cycling and other ecosystem processes, users should use the product's uncertainty with caution and consider inflating the existing LAI product uncertainty when using it within data assimilation frameworks or other studies in which the uncertainty plays an important role.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/12/19/3122/s1,Table S1: Detailed forest community information in the southern China region, Table S2: The detailed information for field LAI measurements at each site, Table S3: Distribution of the ratio of EO uncertainty estimates to bias, based on a comparison against in situ LAI measurements, Figure S1 Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 20

Figure 1 .
Figure 1.Map of the study sites in Southern China.

Figure 1 .
Figure 1.Map of the study sites in Southern China.

Figure 2 .
Figure 2. Time series of the LAI for the field, GEOV3 300 m, GEOV2 1 km, and the MODIS 500 m for different sites.The timespan of the LAI time series is from 2005 to 2017.Each field data point represented the mean value of the measurements from the sample subplots in each month."N" denoted the number of total sample measurements for each site, including the repeat measurements in each month at each subplot.

Figure 2 .
Figure 2. Time series of the LAI for the field, GEOV3 300 m, GEOV2 1 km, and the MODIS 500 m for different sites.The timespan of the LAI time series is from 2005 to 2017.Each field data point represented the mean value of the measurements from the sample subplots in each month."N" denoted the number of total sample measurements for each site, including the repeat measurements in each month at each subplot.

Figure 3 .
Figure 3. Regression between the field LAI and the satellite retrieved LAI.Each point represents the monthly LAI field measurement against the satellite LAI retrieval.The blue solid lines represent the regression between the field LAI against the satellite LAI for all the sites combined on a monthly scale.Point colors and shapes represent different sites.

Figure 3 .
Figure 3. Regression between the field LAI and the satellite retrieved LAI.Each point represents the monthly LAI field measurement against the satellite LAI retrieval.The blue solid lines represent the regression between the field LAI against the satellite LAI for all the sites combined on a monthly scale.Point colors and shapes represent different sites.

Figure 4 .
Figure 4. Distribution of the ratio of EO uncertainty estimates to bias, based on a comparison against in situ LAI measurements.Dashed lines represent the median values and the color represents different EO products.The vertical black line is the reference line where the ratio is equal to 1.A ratio larger than 1 indicates the uncertainty estimate is larger than the bias and is therefore considered robust.

Figure 4 .
Figure 4. Distribution of the ratio of EO uncertainty estimates to bias, based on a comparison against in situ LAI measurements.Dashed lines represent the median values and the color represents different EO products.The vertical black line is the reference line where the ratio is equal to 1.A ratio larger than 1 indicates the uncertainty estimate is larger than the bias and is therefore considered robust.

Figure 5 .
Figure 5. Periodic characteristics of the LAI time series from 2005 to 2017 for different sites.(a): Baseline values of the LAI time series; (b): Phase of the LAI time series; (c): Amplitude of the LAI time series; (d): Relative amplitude of the LAI time series.Point colors and shapes represent different sites.A detailed description of the calculation can be found in the Methods section of this paper.

Figure 5 .
Figure 5. Periodic characteristics of the LAI time series from 2005 to 2017 for different sites.(a): Baseline values of the LAI time series; (b): Phase of the LAI time series; (c): Amplitude of the LAI time series; (d): Relative amplitude of the LAI time series.Point colors and shapes represent different sites.A detailed description of the calculation can be found in the Methods section of this paper.

:
Forest appearance for each forest in this study in the southern China region, Figure S2: TOPO30 DEM (1 km) in the China region.
Figure S3: Location of the sampling subplots and the GEOV2 1km LAI pixels showed on the Landsat EVI scenes for different sites, Figure S4: Location of the sampling subplots and the GEOV2 1 km LAI pixels showed on the SARTM 90 m DEM images for different sites.
Figure S5: Frequency distribution of the ratio of LAI products uncertainty with LAI products' bias against the in situ measurements for different sites.
Figure S6: Regression analyses between the Landsat EVI/NDVI values and the field observed LAI values at sampling subplots, Figure S7: Regression analyses between the proportion of forests cover (LAI pixel level) and the remotely sensed LAI bias against the field measured LAI for different LAI products, Figure S8: Summary of the dates for the MODIS LAI pixels which were influenced by the cloud in this study.

Table 1 .
Background information for each site.MAT: mean annual temperature; AP: annual precipitation.

Table 2 .
Mean site LAI values (m 2