Evaluating Prediction Models for Mapping Canopy Chlorophyll Content Across Biomes

: Accurate measurement of canopy chlorophyll content (CCC) is essential for the understanding of terrestrial ecosystem dynamics through monitoring and evaluating properties such as carbon and water ﬂux, productivity, light use e ﬃ ciency as well as nutritional and environmental stresses. Information on the amount and distribution of CCC helps to assess and report biodiversity indicators related to ecosystem processes and functional aspects. Therefore, measuring CCC continuously and globally from earth observation data is critical to monitor the status of the biosphere. However, generic and robust methods for regional and global mapping of CCC are not well deﬁned. This study aimed at examining the spatiotemporal consistency and scalability of selected methods for CCC mapping across biomes. Four methods (i.e., radiative transfer models (RTMs) inversion using a look-up table (LUT), the biophysical processor approach integrated into the Sentinel application platform (SNAP toolbox), simple ratio vegetation index (SRVI), and partial least square regression (PLSR)) were evaluated. Similarities and di ﬀ erences among CCC products generated by applying the four methods on actual Sentinel-2 data in four biomes (temperate forest, tropical forest, wetland, and Arctic tundra) were examined by computing statistical measures and spatiotemporal consistency pairwise comparisons. Pairwise comparison of CCC predictions by the selected methods demonstrated strong agreement. The highest correlation (R 2 = 0.93, RMSE = 0.4371 g / m 2 ) was obtained between CCC predictions of PROSAIL inversion by LUT and SNAP toolbox approach in a wetland when a single Sentinel-2 image was used. However, when time-series data were used, it was PROSAIL inversion against SRVI (R 2 = 0.88, RMSE = 0.19) that showed greatest similarity to the single date predictions (R 2 = 0.83, RMSE = 0.17 g / m 2 ) in this biome. Generally, the CCC products obtained using the SNAP toolbox approach resulted in a systematic over / under-estimation of CCC. RTMs inversion by LUT (INFORM and PROSAIL) resulted in a non-biased, spatiotemporally consistent prediction of CCC


Introduction
Canopy chlorophyll content (CCC), the sum of chlorophyll a and b in a group of foliage per unit ground area, is a biodiversity variable highly related to ecosystem functioning and helps to understand the dynamics of photosynthesis, plant responses to environmental change, variation in genetics, and diversity in ecology [1]. CCC is an important indicator of ecosystem health and vegetation physiological status [2]. It is an essential variable in monitoring the terrestrial carbon budget by supporting an accurate estimate of primary productivities [3,4] used as an input parameter for biosphere models to predict carbon and water changes [5], and light use efficiency [6]. Information on the amount and distribution of CCC has been utilized to answer many ecological questions related to monitoring and evaluating terrestrial ecosystem properties such as identifying types of vegetation, mapping vegetation cover, and understanding the condition of vegetation [7]. Changes in CCC indicate the effects of disease, nutritional, and environmental stresses [8][9][10][11]. CCC can also be used for forage quality assessment, ecosystem classification, and biomass estimation, as well as being a key input to estimate the indicators for the Convention on Biological Diversity (CBD) such as trends in carbon stocks and patterns in resilience within ecosystems of Aichi target 15, and net primary productivity of Achi target 3 [12]. The spatially and temporally contiguous information on CCC can also be used to support measuring indicators of the sustainable development goals (SDGs) related to the sustainable use of terrestrial ecosystems.
Remote sensing provides the opportunity to drive comprehensive variables for assessing and monitoring biodiversity globally [13] and helps to fill the spatial and temporal gaps left by in situ observations. Thus, remote sensing data provide a viable solution to predict CCC and minimizes the effort by ecologists to collect chlorophyll data through laborious field surveys, cover relatively small extents and short temporal periods to track biodiversity changes.
There is a sturdy correlation between canopy reflectance and chlorophyll that could be used to estimate CCC from imagery. Recent advances in satellite sensors deliver increased resolution and sensitivity of spectral reflectance to CCC. Likewise, a multitude of vigorous methods has been established that correlate CCC and earth observation data [14][15][16][17][18]. The reflectance spectrum between 680 and 760 nm acquired using remote sensing sensors, which is referred to as the red edge position (REP), has been broadly utilized to predict leaf and canopy chlorophyll content [19,20]. This is because a surge in chlorophyll content leads to a shift in the absorption feature towards longer wavelengths [21]. In the literature, the approaches used to predict CCC from earth observation data are categorized into two broad classes: (1) empirical-based; and (2) physical-based [22,23]. Verrelst, et al. [24] further extended the algorithm classes into subclasses and combinations thereof: (i) parametric regression, (ii) non-parametric, (iii) physically-based, and (iv) hybrid methods. Each group of algorithms has its advantages and drawbacks, which are discussed in Verrelst, Camps-Valls [24].
Irrespective of the methods, the prediction accuracy of remotely sensed CCC is affected by factors such as soil type, the existence of non-green components, canopy structure, and shadows [25,26]. Moreover, differences in leaf internal structure, thickness, water content, leaf area index (LAI), foliage clumping, canopy closure, and background in different biomes may alter the relationship between CCC and reflectance [27,28]. Heterogeneity of land surface texture is another source of error when locally developed algorithms up-scaled to regional and global scales [29]. The absence of good-quality high-resolution imagery covering the entire Earth from one sensor may be solved by combining data obtained from multiple sensors, which increases the challenges of global mapping [30]. It is likely that leaves with the same reflectance spectra may have different reflectance factors in different forest canopies [31,32]. As a result transferring predictive models developed for one biome, universally across all biomes in broad spatial extents containing different species or plant functional types is challenging. Previous studies showed that CCC could be retrieved from hyperspectral remote sensing with high accuracy at a relatively small scale, but as soon as the observation level changes to landscape, the accuracy tends to weaken [33]. Algorithms initially developed at a local level are particularly affected by such additional heterogeneity dynamics if extrapolated over a larger extent [25,33,34].
Hence, CCC mapping across biomes for biodiversity monitoring in terrestrial ecosystems requires a reliable operational mapping strategy to be developed. A diversity of mapping schemes and classification methods were proposed for the global mapping of biophysical variables such as LAI and land cover products [35,36]. So far, the only available remote sensing-based chlorophyll products are from ENVISAT MERIS satellite data at 300 m spatial resolution including the MERIS terrestrial chlorophyll index (MTCI) [37], which is later modified to Sentinel-3 Ocean land color imager Terrestrial Chlorophyll Index (OTCI) [38], and the recently developed global distribution of leaf chlorophyll content with canopy and leaf-level radiative transfer modeling procedures [39]. To our knowledge, there are no algorithms validated across biomes for global mapping of CCC from high-resolution satellite imagery. Ali, et al. [40] has recently identified four potential methods that can be used for wide extent retrieval of CCC through a review of the-state-of-art algorithms and validated using in situ data namely simple ratio vegetation index (SRVI) and partial least square regression (PLSR) from statistical-based algorithms, and two physical-based algorithms: INFORM/PROSAIL inversion using a look-up table and the biophysical retrieval tool integrated into the Sentinel application platform (SNAP toolbox), which is a PROSAIL inversion using artificial neural network (ANN) (hereafter referred to as SNAP toolbox approach). To identifying the operationally feasible algorithm(s) that can be upscaled to regional and global levels demands evaluating their accuracy, temporal stability, and transferability across different biomes. Therefore, this study evaluated the robustness and spatiotemporal consistency of those methods proposed by Ali, Darvishzadeh [40] for large scale mapping of CCC for biodiversity monitoring in terrestrial ecosystems. Specifically, we aimed at pairwise comparison of the four methods performance in different pilot sites of major biomes in retrieving CCC from freely available satellite observations, focusing on Sentinel-2.

Pilot Sites
To verify and validate the adequacy of the proposed algorithm(s) for monitoring CCC at a large scale, the selected pilot sites should represent the major biomes of terrestrial ecosystems. Four Biomes (i.e., (1) Boreal Taiga and Artic Tundra, (2) Wetlands, (3) Temperate/Mediterranean forests, and (4) Tropical/sub-tropical rain forest) with diverse vegetation types across gradients of climate, altitude, and latitude were considered. Each biome is represented by one pilot site ( Figure 1). The selected four pilot sites from four terrestrial biomes were:

1.
Kytalyk (Arctic tundra) is located south of the East Siberian Sea coast, Russia (79.82 • N, 147.47 • E). It is part of a low-Arctic tundra nature reserve in the Indigirka lowlands. The site is characterized by very thick (>100-meter depth) ice-rich permafrost [41]. Kytalyk is a habitat of different tundra vegetation types, including tussock sedge, short shrubs, and moss tundra [42].

2.
La Camargue (wetland ecosystem) is one of the Ramsar sites, which is a biosphere reserve in the Rhône delta in Southern France (43.53 • N, 4.50 • E). The major natural habitats of La Camargue contain lagoons, brackish/freshwater marshes, halophilous scrubs, and steppes. In this ecosystem, rice and irrigated crops are predominantly intermingled with the natural wetland vegetation.
Camargue is a species-rich Mediterranean wetland with more than 1200 plant species [43]. The quality and quantity of water that is available year-round highly influence the functional diversity of the pilot site. Significant parts of Camargue naturally dry up during the summer season [44].

Satellite Imagery
To assess and compare the robustness and spatiotemporal consistency of the chosen methods in predicting CCC from remote sensing data for global prediction, Sentinel-2 images of the pilot areas were utilized. Sentinel-2 Level 1C products acquired in 2017 and/or 2018 on a relatively cloud-free

Satellite Imagery
To assess and compare the robustness and spatiotemporal consistency of the chosen methods in predicting CCC from remote sensing data for global prediction, Sentinel-2 images of the pilot areas were utilized. Sentinel-2 Level 1C products acquired in 2017 and/or 2018 on a relatively cloud-free date(s) (< 10% cloud cover) in each pilot site were used. As indicated in Table 1, 24 Sentinel-2 images in total were downloaded from the ESA Copernicus open access hub (https://scihub.copernicus.eu) for this study. Sentinel-2 has three sets of spatial resolution (i.e., 10 m, 20 m, and 60 m). Spectral reflectance in the visible and NIR and SWIR spectral regions are acquired in four, six, and three bands at each spatial resolution, respectively ( Table 2). The Sen2cor vs. 2.5.5 stand-alone software, which is freely distributed under the GNU general public license (http://step.esa.int/main/third-party-plugins-2/ sen2cor/), was used to perform atmospheric correction and produce surface reflectance (level 2A product) from the top-of-atmosphere (TOA) reflectance of Sentinel-2 (level 1c product). The level 2A product was resampled to 20 m spatial resolution, and spectral information in only ten bands (bands 2, 3, 4, 5, 6, 7, 8, 8a, 11, and 12) were used in this study.

Selected Algorithms
We investigated the robustness and spatiotemporal consistency of various approaches to understanding their potential, limitations, and transferability across biomes. More details about the following selected methods can be found in Ali, Darvishzadeh [40].
(a) Simple ratio vegetation indices Two simple ratio vegetation indices optimized for forests and non-forest vegetation were used. The SRVIs use spectral information in two spectral bands to compute CCC. During experimental analysis, Ali, Darvishzadeh [40] revealed that the simple ratio vegetation index 1 (SRVI_1) of Sentinel-2 band 8a (865 nm) and band 4 (665 nm) as a good predictor of CCC in a forest ecosystem (Equation (1)). Another simple ratio vegetation index (SRVI_2) of Inoue, Guérif [10], which is calibrated and validated for retrieving CCC from remote sensing data for a wide range of crops and natural grasses (Equation (2)) was used for non-forest vegetation.
where R865, R835, R704, and R665 are reflectances in the center wavelengths of the Sentinel-2 band setting. CCC was then generated by applying the fitting equations (Equations (3) and (4)) proposed for the two SRVIs on Sentinel-2 level 2A product (TOC reflectance). Among the non-parametric regression methods, the PLSR was utilized to generate CCC products at the four pilot sites. We used the PLSR trained on a spectral subset of eight bands (B2, B3, B4, B5, B6, B8a, B11, and B12) of Sentinel-2 with five components, which was found to be the accurate setup for CCC retrieval from Sentinel-2 data after testing different subsets of spectral information and optimizing the number of components of explanatory variables. The selected bands were those which showed higher variables importance in the projection (≥1) during the PLSR training.
(c) INFORM and PROSAIL radiative transfer models inversion using look-up table (LUT) The CCC products predicted by radiative transfer models (RTMs) inversion were performed by coupling the PROSPECT leaf model with two canopy models: Scattering by Arbitrarily Inclined Leaves (SAIL), and Invertible Forest Reflectance model (INFORM). The two RTMs were used to simulate and generate LUTs for non-forest and forest vegetation pilot sites' Sentinel-2 spectra. Then merit functions applied on the LUTs to retrieve the CCC products. The essential elements of the parametrization, LUT generation, and inversion using the INFORM and PROSAIL models are described briefly below.
(i) Parameterization and generation of LUT using INFORM INFORM [27,48] is parameterized by leaf, canopy, and sensor parameters ( range. It is always a trade-off between the size of the LUT and the computation time of the inversion. The larger the size of the LUT, the higher the chance of the simulated spectra contain all possible combinations of the input parameters, but the inversion becomes computationally expensive. Therefore, a LUT of 200,000 spectra was used, which was a recommendation of several authors [49,50]. To account for model uncertainties and reduce auto-correlation between the spectrum and input variables, a random Gaussian noise value of 0.3% was added to each simulated spectrum. (ii) Parameterization and generation of LUT using PROSAIL PROSAIL is used for the generation of the canopy reflectance spectra of 'short vegetation' such as wetlands, taiga, and tundra. Spectral simulation using PROSAIL requires leaf, canopy, and sensor configuration parameters. A list of the parameters used and their range based on prior knowledge in the literature are presented in Table 4.
A LUT of 100,000 spectra were built by varying the inputs randomly within their range. This size of LUT has been confirmed to be large enough for retrieval of vegetation properties in different vegetation [50][51][52][53]. Similar to the INFORM spectra, a random Gaussian noise value of 0.3% was added to each simulated spectrum.  For both LUT generated by INFORM and PROSAIL, LUT inversion demands searching for each measured spectra (Sentinel-2) it is matching from the simulated spectra (INFORM or PROSAIL). The least root mean square error (RMSE) was the criterion to find the matching spectra through a comparison of the measured and simulated spectra (Equation (5)).
where R measured is a Sentinel-2 reflectance at wavelength λ, and R simulated is a spectrum at wavelength λ in the LUT, and n stands for the number of wavelengths. The LUT based CCC retrieval was performed using three Sentinel-2 bands, which are located in the red, red-edge, and NIR transition zone. This region of the electromagnetic spectrum known for its sensitivity for subtle variation in chlorophyll. The group of input variables that provided the spectra in the LUT with the least RMSE was taken as the best estimates. However, simplification of reality in the RTM and parameterization, as well as calibration and atmospheric correction errors in the remote sensing dataset, the lowest RMSE criterion might not certainly offer the best solution. Consequently, for each Sentinel-2 (measured) spectrum, the 100 similar spectra with relatively smaller RMSE were obtained from the LUT. From the 100 available solutions (q), the median spectrum CCC value was taken as a final solution.
(d) SNAP toolbox approach Baret [56] proposed retrieval of vegetation biophysical variables by training the ANN on PROSAIL simulated spectra and inverting on atmospherically corrected (surface reflectance) data of Sentinel-2. This approach is grouped under the category of 'hybrid' model by Verrelst, Camps-Valls [24], and has been implemented in the Sentinel Application Platform (SNAP) toolbox. The method uses normalized Sentinel-2 data in eight bands (B3 -B7, B8a, B11, and B12), as well as sensor observation and illumination angles, and predicts vegetation biophysical variables such as LAI, fraction of absorbed photosynthetically active radiation, fraction of vegetation cover, CCC, and water content together with a quality flag product for each SNAP biophysical product. The CCC quality flag was used to mask out pixels with input and/or output out of range values. The unit of the CCC product in SNAP is µg/cm 2 ; however, in our study, we applied a conversion factor to convert it to g/m 2 . See Baret [56] for details of the approach.

Assessment of Methods Transferability
The traditional validation of remote sensing methods requires the collection of in situ data concurrent with the acquisition of remote sensing data. Acquiring such in situ dataset from a range of biomes representing a logical subset of the whole terrestrial ecosystems is very expensive and time-consuming, and thus not feasible. To overcome this limitation and quantitatively compare the selected algorithms for regional and global scale CCC mapping, the following verification (validation) strategies were used.
(a) Spatial distribution consistency The spatial distribution of CCC maps produced using the approaches described in Section 2.3.1 was visually and quantitatively investigated for each pilot site in the four biomes. The ranges of the generated CCC products were compared to the expected CCC range in each biome.
A random sample up to 200 pixels was extracted from each CCC product and used in a paired t-test, (assuming the null hypothesis that pairs maps were identical to each other regardless of which algorithm was used) to examine whether the mean difference in CCC from one map to another was different than would be expected by chance alone. We also applied the two-sample Kolmogorov-Smirnov test [59], which is a non-parametric hypothesis test that evaluates the difference between two distributions, to measure disparities among pairs of CCC products.
(b) A measure of agreement among pairwise CCC products The closeness of the CCC predicted values by pairs of methods was evaluated by computing the strength of the correlation, prediction errors, and precision. The randomly extracted 200 samples of CCC values from each map were used to calculate the coefficient of determination (R 2 ) (Equation (6)), root mean square error (RMSE) (Equation (7)), and Bias (Equation (8)). The more similar pairs of products are the ones with higher R 2 , lower RMSE, and bias close to zero.
where y and y' are the predicted values of the two methods for sample i, and n is the number of samples considered.

(c) Temporal consistency
To examine the temporal consistency and robustness of the candidate algorithms, predictions of CCC were performed using time series Sentinel-2 data. It was evaluated whether the closeness of the predicted CCC values by a pair of methods significantly change through time. Therefore, the candidate algorithms were applied on cloud-free time-series Sentinel-2 data available for the period June 2017 to September 2018 for two of the pilot sites (BFNP and La Camargue) assuming the temporal consistency in these two sites represent the performance of the methods across tall (forest) and short (shrub and grassland) biomes. Sentinel-2 images acquired on the seven dates were used for BFNP. Whereas for La Camargue, fifteen dates of Sentinel-2 images were used. The statistical measures such as R 2 , RMSE, and bias (Equations (6)-(8)) were computed by using several dates of CCC values predicted by a pair of methods to assess how the relationship between pairs of algorithms change over time.

Spatial Distribution Consistency
The expected and predicted CCC ranges are presented in Table 5. The INFORM and PROSAIL ranges are closer to expectations in three of the four biomes (temperate, wetland, and tundra). But, in tropical rain forests, it is the range obtained by the SNAP toolbox which is closest to expectations. Figures 2-5 details the spatial distribution consistency of the CCC predictions by the candidate approaches in the four biomes. Visually, the CCC maps generated by all methods showed a similar pattern in the BFNP (Figure 2). In this pilot site, generally, deciduous forests (mostly beech) showed higher CCC values. However, there are differences among the CCC products of the selected methods for the other three biomes. It is apparent from closer inspection of the four CCC maps generated for the tropical forest ( Figure 3) that there is a more significant variation among the produced products compared to BFNP. The prediction made by the SNAP toolbox indicated higher CCC values in most parts of the Lambir pilot site (tropical forest) than other approaches. The CCC product from the PLSR method deviates from the other three approaches when applied in the Mediterranean wetland ecosystem ( Figure 4) and tundra pilot sites ( Figure 5), and thus no further statistical test performed for PLSR.
Remote Sens. 2020, 12, 1788 20 of 32 prediction made by the SNAP toolbox indicated higher CCC values in most parts of the Lambir pilot site (tropical forest) than other approaches. The CCC product from the PLSR method deviates from the other three approaches when applied in the Mediterranean wetland ecosystem ( Figure 4) and tundra pilot sites (Figure 5), and thus no further statistical test performed for PLSR.    The two-sample Kolmogorov-Smirnov test among pairs of CCC products from INFORM, SRVI, and SNAP did not show significant distribution disparities in the temperate forest pilot site (Table 5). The Violin plot in Figure 6 depicts the similarities and differences of CCC products generated by the selected methods. In contrast, all of the pairs of CCC products showed distribution disparity in tropical rainforest, as shown in Table 6 and the violin plots in Figure 6c. Pairwise comparison of CCC predicted by PROSAIL and SRVI in La Camargue (wetland) exhibited similar distribution. In The two-sample Kolmogorov-Smirnov test among pairs of CCC products from INFORM, SRVI, and SNAP did not show significant distribution disparities in the temperate forest pilot site (Table 5). The Violin plot in Figure 6 depicts the similarities and differences of CCC products generated by the selected methods. In contrast, all of the pairs of CCC products showed distribution disparity in tropical rainforest, as shown in Table 6 and the violin plots in Figure 6c. Pairwise comparison of CCC predicted by PROSAIL and SRVI in La Camargue (wetland) exhibited similar distribution. In contrast, non-significant distribution disparity was observed between CCC products from the SNAP toolbox and SRVI in tundra biome (Table 5). However, the latter two methods predicted CCC values ranges in tundra biome were far beyond the expected range in Table 5.
Remote Sens. 2020, 12, 1788 26 of 32 consistency of the robustness of the methods through time. The seasonal variability of CCC in broadleaf and conifer stands of the Bavarian temperate forest is demonstrated in Figure 8. As expected, the variation in CCC is minimal in conifer stands compared to broadleaf.

The Agreement of CCC Values Predicted by the Selected Methods
We verified how the CCC products generated by the methods chosen are close to each other across the four biomes. As illustrated in Figure 7, the CCC predicted by SRVI agrees with radiative transfer models (RTMs) inversion (INFORM inversion in forests or PROSAIL in short vegetation). The SRVI and RTMs inversion pairs had a higher correlation, lower RMSE, and the bias was close to zero in all the four pilot sites. However, the highest correlation (R 2 = 0.93) was recorded for the SNAP toolbox approach with PROSAIL inversion (Figure 7i).

Temporal Consistency
The selected methods were applied to a time series of Sentinel-2 imagery for the Bavarian Forest pilot site, as well as the short vegetation biome (La Camargue pilot site) in order to assess the The predictions made using the SNAP toolbox approach [56] showed a tendency of overestimation, particularly in a wetland, when compared to SRVI and RTM inversion by LUT. There are weaker agreements between predictions in tundra than any other biome (Figure 7j).

Temporal Consistency
The selected methods were applied to a time series of Sentinel-2 imagery for the Bavarian Forest pilot site, as well as the short vegetation biome (La Camargue pilot site) in order to assess the consistency of the robustness of the methods through time. The seasonal variability of CCC in broadleaf and conifer stands of the Bavarian temperate forest is demonstrated in Figure 8. As expected, the variation in CCC is minimal in conifer stands compared to broadleaf.
consistency of the robustness of the methods through time. The seasonal variability of CCC in broadleaf and conifer stands of the Bavarian temperate forest is demonstrated in Figure 8. As expected, the variation in CCC is minimal in conifer stands compared to broadleaf.  Plotting the predictions made by the candidate algorithms in BFNP against each other for all the available image dates (seven dates) for randomly selected pixels to check if the relationship changes through time did not show a significant difference (Figure 9). The time series scatter plots have a similar pattern to the result obtained by applying the methods in single time Sentinel-2 data demonstrated in Figure 7a-c. Nonetheless, the INFORM inversion against the SNAP toolbox approach (Figure 9c) elucidated a higher correlation to the single date product (Figure 7c) than others (Figure 9a,b).
available image dates (seven dates) for randomly selected pixels to check if the relationship changes through time did not show a significant difference (Figure 9). The time series scatter plots have a similar pattern to the result obtained by applying the methods in single time Sentinel-2 data demonstrated in Figure 7a-c. Nonetheless, the INFORM inversion against the SNAP toolbox approach (Figure 9c) elucidated a higher correlation to the single date product (Figure 7c) than others (Figure 9a,b).   Plotting the predictions made by the candidate algorithms in BFNP against each other for all the available image dates (seven dates) for randomly selected pixels to check if the relationship changes through time did not show a significant difference (Figure 9). The time series scatter plots have a similar pattern to the result obtained by applying the methods in single time Sentinel-2 data demonstrated in Figure 7a-c. Nonetheless, the INFORM inversion against the SNAP toolbox approach (Figure 9c) elucidated a higher correlation to the single date product (Figure 7c) than others (Figure 9a,b).

Discussion
Very little was found in the literature on the global canopy chlorophyll content (CCC) products from high spatial resolution remote sensing data. The present study was designed to determine the best-practice approach that may be implemented for mapping such CCC products through spatiotemporal consistency and robustness comparison across biomes.
One interesting finding is the spatial distribution of estimated CCC by the candidate algorithms displayed similar patterns in the BFNP, where calibration and validation of the methods were performed using an in situ dataset ( Figure 2). Likewise, the violin plot ( Figure 6) and the two-sample Kolmogorov-Smirnov test ( Table 6) among pairs of CCC products proved the absence of significant spatial distribution disparities in this pilot site. These results are in line with previous findings [23,29], which showed the vigor of models when they are tuned for a particular biome.
However, some of the methods, mainly, PLSR, poorly performed in other biomes (Figures 4  and 5). The possible explanation for the weak performance of the PLSR could be due to overfitting of non-parametric methods that stem from excessively multifaceted models that seize more than the underlying correlation [24,64]. This suggests that non-parametric methods are not well suited for global retrieval schemes. Although the PLSR performed well during validation using in situ data in Ali, Darvishzadeh [40], it posed a question about the transferability of the method across biomes. It demands calibration of PLSR for different biomes separately, and this is not feasible if it has to be implemented on an operational basis. Therefore, PLSR was not found suitable for large scale mapping, and thus we compared to the other three algorithms (i.e., the SNAP toolbox, SRVI, and INFORM/PROSAIL) for their operational feasibility through observing spatiotemporal consistency and robustness across biomes. Consistent with the literature, this research found that statistical methods are site and sensor-specific and rely on sampling circumstances, and may vary in space and over time [29,65].
The predictions made using SRVI and the SNAP toolbox approach, which is a combination of statistical and physical models, resulted in a systematic over/under-estimation of CCC when applied in different biomes. For instance, the SNAP toolbox approach [56] showed a tendency of overestimation, particularly in the wetland study site, when compared to SRVI and RTM inversion by LUT (Figure 7g-i). This result largely supports the finding of other authors [66][67][68] who highlighted the limitation of statistical methods in predicting vegetation variables from earth observation data.
In terms of predictive values agreement, all pairs of CCC products produced by the candidate approaches showed a good correlation. The highest correlation (R 2 = 0.93) was observed between predictions of the SNAP toolbox and PROSAIL inversion by LUT (Figure 7i), and the weakest relationship (R 2 = 0.57) was recorded in the tundra (Figure 7j). The weak agreements between predictions in the tundra may be partly due to the low chlorophyll content of tundra biomes, which is close to the minimum detectable CCC by several of the candidate methods. Weak relationships (R 2 ≤ 0.43) between pigment content and spectral indices have also been reported by Beamish, et al. [69] in a low arctic tundra terrestrial ecosystem.
CCC predictions by INFORM and PROSAIL inversion by LUT exhibit ranges closer to expectations in three of the four biomes (Table 5), which implies that RTM approaches are rigorous. In many cases, lower disparities, higher R 2 , lower RMSE, and bias close to zero were observed when other methods compared with RTM inversions by LUT (Table 6, and Figure 7c,e,h,l), and reaffirmed the RTM based approaches applicability in different biomes. It is encouraging to compare this figure with that found by Ali, Darvishzadeh [40] who found that unlike the SRVI and the SNAP toolbox approach, RTM inversions by LUT provided a non-biased prediction.
Temporal consistency verification in the two pilot sites (Bavaria and La Camargue) also portrayed the robustness of the RTM based approaches. Pairwise comparison of the time series analysis products from INFORM inversion against SNAP toolbox approach (R 2 = 0.88, RMSE = 0.25 g/m 2 ) in temperate forest (Figure 9c), and PROSAIL inversion against SRVI (R 2 = 0.88, RMSE = 0.19 g/m 2 ) in the wetland site (Figure 10b) elucidated more similarity to the single date predictions shown in Figure 7c (R 2 = 0.82, RMSE = 0.3 g/m 2 ) and Figure 7h (R 2 = 0.83, RMSE = 0.17 g/m 2 ) respectively than others pairs, which reaffirms the fact that RTM based predictions are spatiotemporally consistent. These results are in accord with a recent study indicating that an RTM inversion is the robust method in the global mapping of the distribution of leaf chlorophyll content from ENVISAT MERIS full resolution (300 m) satellite imagery [39].
However, RTM inversion using LUTs requires a large set of input parameters at the leaf and canopy level to parameterize the models, and the models can be computationally expensive to execute. For instance, the inversion of INFORM took two solid days (result not shown) per pilot site. In contrast, the approach already implemented in the ESA SNAP toolbox and the parametric regression approach (i.e., SRVI) is easy and convenient to retrieve CCC from Earth observation data. It took less than one minute per pilot site. Thus, the SRVI and the SNAP toolbox approach may be the operationally more feasible approaches in terms of computation for the prediction of CCC from large remote sensing datasets. There is, therefore, a definite need for a trade-off between the accurate prediction of CCC globally and computational efficiency. Recent advances in computer science and cloud computing systems tremendously improved computational efficiency. This implies that the computationally demanding nature of RTMs inversion does not preclude them from being operationally feasible algorithms.

Conclusions
The purpose of this work was to determine the best-practice approach that can be implemented to map canopy chlorophyll content (CCC) across biomes from high spatial resolution remote sensing data. This study builds upon our earlier study, where the selected methods have been validated using existing in situ data [40]. In this study, we further compared their spatiotemporal consistency and robustness across different biomes to recommend the method(s) that suit for large-scale mapping of CCC from remote sensing data.
The major finding was that all methods exhibit spatiotemporal consistency. However, the statistical (e.g., SRVI) and hybrid (e.g., SNAP toolbox) approaches lead to a systematic over/under-estimation of CCC when applied in different biomes. CCC predictions by INFORM and PROSAIL inversion by LUT were rigorous and much closer to expected ranges. They were found non-biased and robust predictor across biomes. Therefore, based on the theoretical analysis of the rewards and drawbacks of the algorithms, spatiotemporal consistency across biomes, and robustness comparisons, the RTM inversion using LUT approach particularly INFORM for 'forest' and PROSAIL for 'short vegetation' ecosystems are recommended for large scale mapping of CCC from Sentinel-2 data. Future endeavors may include more rigorous validation of the CCC products obtained by the recommended approaches using field data collected for different terrestrial biomes.