Retrieval of the Fraction of Radiation Absorbed by Photosynthetic Components (FAPARgreen) for Forest Using a Triple-Source Leaf-Wood-Soil Layer Approach

The fraction of absorbed photosynthetically active radiation (FAPAR) is generally divided into the fraction of radiation absorbed by the photosynthetic components (FAPARgreen) and the fraction of radiation absorbed by the non-photosynthetic components (FAPARwoody) of the vegetation. However, most global FAPAR datasets do not take account of the woody components when considering the canopy radiation transfer. The objective of this study was to develop a generic algorithm for partitioning FAPARcanopy into FAPARgreen and FAPARwoody based on a triple-source leaf-wood-soil layer (TriLay) approach. The LargE-Scale remote sensing data and image simulation framework (LESS) model was used to validate the TriLay approach. The results showed that the TriLay FAPARgreen had higher retrieval accuracy, as well as a significantly lower bias (R2 = 0.937, Root Mean Square Error (RMSE) = 0.064, and bias = −6.02% for black-sky conditions; R2 = 0.997, RMSE = 0.025 and bias = −4.04% for white-sky conditions) compared to the traditional linear method (R2 = 0.979, RMSE = 0.114, and bias = −18.04% for black-sky conditions; R2 = 0.996, RMSE = 0.106 and bias = −16.93% for white-sky conditions). For FAPAR that did not take account of woody components (FAPARnoWAI), the corresponding results were R2 = 0.920, RMSE = 0.071, and bias = −7.14% for black-sky conditions, and R2 = 0.999, RMSE = 0.043, and bias = −6.41% for white-sky conditions. Finally, the dynamic FAPARgreen, FAPARwoody, FAPARcanopy and FAPARnoWAI products for a North America region were generated at a resolution of 500 m for every eight days in 2017. A comparison of the results for FAPARgreen against those for FAPARnoWAI and FAPARcanopy showed that the discrepancy between FAPARgreen and other FAPAR products for forest vegetation types could not be ignored. For deciduous needleleaf forest, in particular, the black-sky FAPARgreen was found to contribute only about 23.86% and 35.75% of FAPARcanopy at the beginning and end of the year (from January to March and October to December, JFM and OND), and 75.02% at the peak growth stage (from July to September, JAS); the black-sky FAPARnoWAI was found to be overestimated by 38.30% and 28.46% during the early (JFM) and late (OND) part of the year, respectively. Therefore, the TriLay approach performed well in separating FAPARgreen from FAPARcanopy, which is of great importance for a better understanding of the energy exchange within the canopy. Remote Sens. 2019, 11, 2471; doi:10.3390/rs11212471 www.mdpi.com/journal/remotesensing Remote Sens. 2019, 11, 2471 2 of 19


Introduction
The fraction of absorbed photosynthetically active radiation (FAPAR) is a significant biochemical and physiological variable used in tracing the exchanges of energy, mass, and momentum, and is also widely used in many climate, ecological, biogeochemical, agricultural, and hydrology models [1,2].FAPAR is, therefore, an important input parameter and widely used in satellite-based Production Efficiency Models (PEMs) [3][4][5][6] to estimate gross primary productivity (GPP) or net primary production (NPP).
In general, the FAPAR inversion algorithms could be divided into two types: empirical statistical models based on vegetation indexes and physical methods based on the canopy radiation transfer model.Although the empirical statistical model based on vegetation indexes is relatively simple, involves only a few parameters and has high computational efficiency, it is subject to many uncertainties due to factors such as the atmospheric environment, vegetation type, and quality of remote sensing data.The physically based methods could be further divided into two categories.The first type is the direct inversion method, which uses the canopy radiation transfer model to link FAPAR with the canopy spectra [7][8][9].For example, the Moderate Resolution Imaging Spectroradiometer (MODIS) algorithm uses the three-dimensional radiation transmission model to invert FAPAR from the bi-directional reflectance [10][11][12].The Joint Research Centre (JRC) FAPAR algorithm is also based on a physical model that uses a continuous vegetation canopy model [13] to link land surface reflectance with FAPAR.However, these methods are mostly based on the radiative transfer model; thus, the inversion process is complicated for retrieval of FAPAR.The main problem with such methods is that it is difficult to overcome the uncertainty caused by model coupling and spatial heterogeneity.The second type of physically based method is the forward modeling method [14][15][16][17][18][19].Most models of this type are based on the gap fraction model, which determines FAPAR according to canopy structure parameters such as LAI and the clumping index.The disadvantage of this approach is that it relies too much on the accuracy of the canopy structure parameters.Furthermore, it is difficult to accurately determine the soil albedo and extinction coefficient, which are also important parameters needed to determine the contribution of multiple scattering between the soil and canopy to FAPAR [17].
The forest vegetation ecosystem plays an important role in the global ecosystem.However, quantifying the temporal variation in FAPAR green for a forest ecosystem represents an important challenge for remote sensing and ecology researchers as it is extremely difficult to measure FAPAR green at large scales over plant growing seasons directly.Also, previous studies have shown that the contribution of woody components is relatively large: for instance, Asner et al. [37] found that stems increased FAPAR canopy by 10-40%.Therefore, the partitioning of absorbed radiation into photosynthetic and non-photosynthetic parts is very important for better modeling of vegetation photosynthesis and energy exchange within the canopy.
Already, some studies have looked at the estimation of FAPAR green from remote sensing data.Hall et al. [38] estimated FAPAR green using a simple linear relationship between FAPAR canopy and LAI green /LAI total (LAI total denotes the total leaf area index including green and senescent leaves, while LAI green represents the green leaf area index).However, this simple partitioning is problematic because the green and woody components within the canopy do not constitute a simple linear mix in terms of radiation transfer.Zhang et al. [39] first retrieved the biophysical and biochemical variables using the modified PROSPECT model coupled with the SAIL-2 model (hereafter called PROSAIL-2 model), and then calculated FAPAR green and FAPAR canopy using the forward simulation approach.However, FAPAR green retrieval using the PROSAIL-2 model is relatively complex and needs several physiological and biochemical parameters as model inputs.Gitelson et al. [40] also separated FAPAR canopy into photosynthetically active green components (FAPAR green ) and non-photosynthetic active components using the ratio LAI green /LAI total for maize and soybeans.The relationship between vegetation indices and FAPAR green was also used to retrieve FAPAR green [41,42].Nevertheless, to date, the current FAPAR green products do not take into account the effect of non-photosynthetic components on canopy radiative transfer.
In this study, we aim to develop an operational algorithm for partitioning FAPAR canopy into FAPAR green and FAPAR woody for forest types.A simple triple-source leaf-wood-soil layer model (TriLay) that describes the radiation transfer within the canopy-soil system is presented.FAPAR canopy is first separated into the fraction of PAR absorbed by the canopy for downwelling radiation (FAPAR canopy↓ ) and the fraction of PAR absorbed by the canopy for the upwelling radiation reflected by the soil background (FAPAR canopy↑ ).Then, FAPAR canopy↓ and FAPAR canopy↑ are further split into the fraction of radiation absorbed by photosynthetic components (FAPAR green ) and that absorbed by non-photosynthetic components (FAPAR woody ) using the TriLay model.Finally, the FAPAR green , FAPAR woody , and FAPAR canopy products are generated using the MODIS albedo (MCD43A3), LAI (MCD15A2H), land cover (MCD12Q1), clumping index (CI), and soil albedo products based on the TriLay approach, and the discrepancies between different FAPAR products are used to investigate the contributions of woody components to the canopy-absorbed radiation.The partitioning of absorbed radiation into green and woody parts using the TriLay model is done not just to provide FAPAR green and FAPAR woody -which is of great importance for better understanding the energy exchange within the canopy.The consideration of woody components should also improve the accuracy of FAPAR green estimates, which is important for better modeling of vegetation photosynthesis.

Satellite Datasets
In order to produce FAPAR green products, several satellite datasets were utilized in this study, including MODIS LAI products, land cover products, CI products, and soil albedo products from 2017.Data simulated by the LESS model [43] was used to validate the retrieved FAPAR green products.A mid-latitude region (tile h10v05, covering 30.00 • N-40.00 • N and 80.00 • W-104.43 • W) was selected to investigate the discrepancy between FAPAR green and other FAPAR products because there were abundant forest vegetation types within this MODIS tile.

MODIS LAI/FAPAR Products (MCD15A2H)
MCD15A2H V006 is a MODIS eight-day composite LAI/FAPAR product that includes FAPAR, LAI, and quality control (QC) data with a resolution of 500 m [44].The main retrieval algorithm for LAI and FAPAR contains a Look-up-Table (LUT) based on a 3D radiation transfer model [21] that uses the atmospherically corrected Red and near-infrared (NIR) Bidirectional Reflectance Function (BRF) [45].A back-up algorithm based on the empirical relationships between the Normalized Difference Vegetation Index (NDVI) and LAI and FAPAR at the canopy scale is used at the same time.Also, for the biome types and typical conditions that are considered, observed, and modeled spectral BRFs and soil patterns are compared for each pixel.The LAI and FAPAR values that lie within a fixed level of uncertainty are then taken to be acceptable.Finally, averaged values of LAI and FAPAR are used as the eventually retrieved values [20].

MODIS Land Cover Product (MCD12Q1)
The MODIS Land Cover Type Product (MCD12Q1) provides land cover maps with a temporal resolution of one year and a spatial resolution of 500m at a global scale from 2001 until the present; it includes several classification schemes (International Geosphere-Biosphere Programme (IGBP), University of Maryland (UMD), LAI, BIOME-Biogeochemical Cycles (DBC), Plant Functional Types (PFT), FAO-Land Cover Classification System land cover (LCCS1), etc.).The main algorithm used in MCD12Q1 is a supervised classification method (decision tree) combined with a boosting technique [20] based on MODIS reflectance data [46,47].The International Geosphere-Biosphere Program (IGBP) classification scheme, which defines 17 land cover types, was used in this study.

Global Clumping Index (CI) Product
The clumping index (CI) signifies the characteristics of groups of foliage in the canopy, and thus, it is a significant parameter describing the structure of the vegetation canopy [48].According to Jiao et al. [49], a method based on the MODIS Bidirectional Reflectance Distribution Function (BRDF) and a linear relationship between the CI and the normalized difference between the angular indexes of the hotspot and dark spot (NDHD) could be utilized to generate CIs within a valid range (0.33 to 1.00); a back-up algorithm is also used to substitute values for the invalid CIs [49].A global CI dataset with an eight-day temporal resolution and 500 m spatial resolution covering the period from 2002 to the present has been produced by Jiao et al. [49].

Global Soil Albedo Product
A non-linear spectral mixture model (NSM) model proposed by Liu & Zhang (2018) [15,50] was used to retrieve the global visible (VIS) soil albedo.The main idea in the NSM model is the dual-source vegetation-soil layer approach.In this approach, it is assumed that the leaves are located in the upper canopy, while the soil components are found in the lower canopy.Based on this assumption, the canopy albedo can be approximated as a non-linear mixture of the "pure" vegetation and soil parts [15].For pixels with abnormal values (smaller than 0.02 or greater than 0.3), prior values acquired by a global database of land surface parameters at 1 km resolution (ECOCLIMAP) [51,52] and the yearly composite value were used instead.Finally, the estimated global VIS soil albedo based on the NSM model is obtained, having a spatial resolution of 500 m and a temporal resolution of eight days.

Data Simulated by the LESS Model
To quantitatively evaluate the performance of the proposed TriLay approach for estimating FAPAR green , the LargE-Scale remote sensing data and image simulation (LESS) framework model [43] was employed to generate a simulated dataset that covered most of the conditions found in forests.LESS is a ray-tracing based 3D radiative transfer model which can simulate remote sensing data and images over large-scale and realistic 3D scenes (http://lessrt.org/).LESS employs a weighted forward photon tracing (FPT) method to simulate multispectral bidirectional reflectance factor (BRF) or flux-related data (e.g., downwelling radiation) and a backward path tracing (BPT) method to generate sensor images (e.g., fisheye images) or large-scale (e.g., 1 km 2 ) spectral images.The accuracy of LESS is evaluated with other models as well as field measurements in terms of directional BRFs and pixel-wise simulated image comparisons, which shows very good agreement [43,53].
The modeling area is located in the Genhe Forestry Reserve (Genhe) (120 • 12 to 122 • 55 E, 50 • 20 to 52 • 30 N), Greater Khingan of Inner Mongolia, Northeastern China.It has a hilly terrain with 75% forest cover, which is mainly composed of Dahurian Larch (Larix gmelinii) and White Birch (Betula platyphylla Suk.).A pure plot of Dahurian Larch (L9) is established and is selected for modeling, with its location in Figure 1.The position, crown width, breast diameter, tree height, and transmittance for trees were entered into the LESS model.A total of eight scenes were constructed, with each individual tree in a scene consisting of branches and leaves.The branches and size of the woody area were kept the same in all the scenes, but the leaf area was changed to produce scenes with different values of the LAI.For each scene, different FAPAR values (including FAPAR canopy , FAPAR green and FAPAR woody ) were calculated under black-sky (the ratio of diffuse light is zero, and nine different solar zenith angles of 0 • -80 • at 10 • intervals were set) and white-sky conditions (the ratio of diffuse light is 1).The main input parameters used in the LargE-Scale remote sensing data and image simulation (LESS) model are listed as Table 1.Finally, a total of 72 black-sky simulations and eight white-sky simulations were achieved for the different conditions giving.

Algorithms for Estimating Global FAPAR green and FAPAR woody Datasets
To split the fraction of radiation absorbed by photosynthetic components (FAPAR green ) from FAPAR canopy , a novel triple-source leaf-wood-soil layer model was proposed for generating global FAPAR green products for forest vegetation types.Figure 2 is a flowchart of the process used to retrieve global FAPAR green products.
First, FAPAR canopy is split into two parts: the fraction of PAR absorbed by the canopy for the downwelling radiation (FAPAR canopy↓ ) and that absorbed by the canopy for the upwelling radiation reflected by the soil background (FAPAR canopy↑ ).Then, FAPAR canopy↓ and FAPAR canopy↑ are further split into the fraction of radiation absorbed by photosynthetic components (FAPAR green ) and that absorbed by non-photosynthetic components (e.g., branches and stems, hereafter called FAPAR woody ).Finally, FAPAR green and FAPAR woody can be calculated separately using the TriLay approach.

The Triple-Source Leaf-Wood-Soil Layer Model
A triple-source leaf-wood-soil layer model (TriLay) was developed to model the radiation transfer within the vegetation-soil system-this model is illustrated as Figure 3.In this study, we used the layer approach to illustrate the distribution of leaves and soil in the whole canopy, which is consistent with the approach used in our previous study [15].The layer approach assumes that the canopy consists only of green components and woody components and that all leaves are found above the canopy.The NSM model [15] was then used to simulate the canopy albedo.
The fraction of PAR absorbed by the canopy, FAPAR canopy , can be separated into the fraction of PAR absorbed by the canopy for the downwelling radiation (FAPAR canopy↓ ) and that for the upwelling radiation reflected from the soil background (FAPAR canopy↑ ): where FAPAR canopy↓ describes the fraction of downwelling radiation absorbed within the canopy assuming the soil background is dark, and FAPAR canopy↑ describes the fraction of upwelling radiation absorbed within the canopy due to the interaction between the ground (soil and understory) and the canopy.A point worth emphasizing is that, assuming a black soil background, the FAPAR canopy↓ is equal to the FAPAR canopy .Therefore, FAPAR canopy↓ can be given by: where τ PAI is the transmittance of the whole canopy, which contains green and woody parts, FVC is the fraction of vegetation cover, and Albedo pure is the visible (VIS) albedo for pure vegetation, which is represented by the VIS albedo for vegetation with a "saturated" LAI value (e.g., LAI ≥ 6) [14].
According to the statistical results by Liu et al.
[15], Albedo pure was set to 0.025 (for white-sky conditions) and 0.020 (for black-sky conditions) for all woody vegetation types (including evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, and deciduous broadleaf forest).τ PAI can be calculated as the product of τ LAI and τ WAI , and the directional transmittance can be determined using the gap fraction model [16,53,54]: where LAI and WAI are the leaf area index and woody area index, respectively, k 1 and k 2 are extinction coefficients for the green components and woody components, respectively, and θ is the solar zenith angle.k 1 can be determined using the leaf absorptance in the VIS band and was set to 0.88 based on simulations made using the PROSPECT-5 model and also the Leaf optical properties experiment 93 (LOPEX'93) and Leaf Optical Properties Database (an experiment conducted at the National Institute for Agricultural Research in Angers, France in June 2003) [55][56][57].The woody components were assumed to be opaque with a constant extinction coefficient (k 2 ) of 0.91 based on the simulated LESS data.G(θ) is the projection of the unit foliage area on the plane perpendicular to the solar incident direction, θ.For green leaves, G(θ) is normally given a value of 0.5 for canopies with a spherical leaf angle distribution.For woody components, we assuming the woody components have the same angular distribution as green leaves, which also assumed by Chen et al. [58], Kucharik et al. [59], and Sea et al. [60].CI is the clumping index; we also assume the same CI values for both green leaves and woody components based on the findings of Chen et al. [61] and Zou et al. [62].
The gap fraction model is also used to calculate FVC for a fixed solar zenith angle of 0 • and a fixed value of 1 for the leaf extinction coefficient: Similarly, FAPAR canopy↑ can be calculated as follows: where Albedo soil is the soil visible albedo, which can be generated by the NSM model proposed by Liu et al. [15].τ ws PAI is the canopy transmittance under white-sky conditions, which can be calculated as the product of the white-sky transmittance of leaves (τ ws LAI ) and the white-sky transmittance of woody components (τ ws WAI ): τ ws PAI = τ ws LAI × τ ws WAI (10)

Determination of Woody Area Index
In general, it is expensive and time-consuming to make accurate estimates of the woody area index (WAI), and destructive sampling is often the only option available for the quantification of the WAI in tropical evergreen forests [63].Therefore, for generating global FAPAR green datasets, the use of accurate estimates of WAI is unrealistic.Hence, in this study, we aimed to determine the global WAI using MCD15A2H LAI data and assumed that the WAI was constant within a given year.
First, we assumed that the woody-to-total area ratio is measured during the peak growth stage (July-August-September, JAS) when the LAI has its maximum value for the year.The woody-tototal area ratio was determined according to the forest types listed in the MODIS land cover product (MCD12Q1), including evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), deciduous needleleaf forest (DNF) and mixed forest (MF).Therefore, WAI values for different forest vegetation types were calculated using a simple linear relationship between the plant area index (PAI) and the WAI: WAI = PAI × ratio woody = LAI max 1−ratio woody × ratio woody (11) where ratio woody is the mean value of the woody-to-total area ratio for various forest types, as given in the literature [58,64,65].LAI max is the maximum value of LAI within a given year; this was acquired from MCD15A2H LAI products.

Separating FAPAR green and FAPAR woody from FAPAR canopy
In a similar way to Equation (1), the fraction of PAR absorbed by the green and woody components can also be separated into the fraction of PAR absorbed by the canopy for the downwelling radiation and that for the upwelling radiation reflected from the soil background.Firstly FAPAR canopy↓ is split into FAPAR green↓ and FAPAR woody↓ : where FAPAR green↓ and FAPAR woody↓ can be obtained as FAPAR woody↓ = FAPAR canopy↓ × ratio woody × w 2↓ (14) where ratio green and ratio woody are the ratio of leaf area index to plant area index and woody area index to plant area index, respectively.w 1↓ and w 2↓ are the weighting coefficients for the green (i.e., photosynthetic) and woody components in terms of the radiation transfer within the canopy.The terms involved in equations ( 10) and ( 11) can be given as: w 1↓ × ratio green + w 2↓ × ratio woody = 1 FAPAR green↓ and FAPAR woody↓ can then be obtained by solving equations ( 13)-( 18 Finally, FAPAR green and FAPAR woody can be calculated as: FAPAR woody = FAPAR woody↑ + FAPAR woody↓ (24)

Validation of the TriLay Method using Simulations made by the LESS Model
It is too challenging to obtain in-situ measurements of FAPAR green for forests, and so the simulated dataset (Table 1) derived using the LESS model was used for the validation of the TriLay model.Seventy-two black-sky simulations of FAPAR canopy , FAPAR green , and FAPAR woody together with eight white-sky simulations were available for validation.Figure 4 illustrates the validation results for FAPAR canopy , FAPAR green , and FAPAR woody .The estimated and simulated FAPAR values are distributed close to the 1:1 line.Also, it can be seen that the TriLay approach can produce accurate estimates of FAPAR canopy , giving Root Mean Square Error (RMSEs) of 0.048 and 0.024 for black-sky and white-sky conditions, respectively, as against the LESS-simulated values.The corresponding R 2 values are 0.945 and 0.999.For FAPAR green , the validation results give RMSEs of 0.064 and 0.025, respectively, for black-sky and white-sky FAPAR.Finally, it can be seen that FAPAR woody can also be accurately estimated: the RMSE and R 2 values are 0.042 and 0.709, respectively, for black-sky conditions, and 0.014 and 0.992 for white-sky conditions.These results show that the TriLay approach can be used to accurately estimate FAPAR canopy , FAPAR green , and FAPAR woody for forest land cover types.

Comparison of Different Methods using the LESS Simulations
Hall et al. [38] estimated FAPAR green based on a simple linear relationship between FAPAR canopy and FAPAR green .They used the ratio LAI green /LAI total to determine FAPAR green : FAPAR green = FAPAR canopy × LAI green total LAI (25) where total LAI is the PAI mentioned above.In order to test the accuracy of the linear mixture method, we also validated the FAPAR green estimated by the linear mixture method [38] using LESS-simulated data.
Figure 5 shows the accuracy assessment results for the linear mixture method.The results show a noticeable underestimation for FAPAR green and a huge overestimation for FAPAR woody .Although the FAPARs retrieved using the linear mixture method are highly correlated with the LESS simulation, with R 2 values of 0.979, 0.996 for FAPAR green , and 0.934, 0.985 for FAPAR woody under black-sky and white-sky conditions, the corresponding RMSE values are much higher than those found using our TriLay approach-0.114,0.106 for FAPAR green as against 0.064, 0.025, and 0.113, 0.106 for FAPAR woody as against 0.042, 0.014 under the black-sky and white-sky conditions, respectively.The FAPAR values derived without considering the woody components (FAPAR noWAI ) were also validated using the LESS simulations; the results are shown in Figure 6.These results give an RMSE of 0.071 and R 2 of 0.920 for the black-sky conditions and the corresponding values of 0.043 and 0.999 for the white-sky conditions.Although the R 2 values are higher, the RMSEs are still greater than those found using the TriLay approach.In order to further compare the performances of different methods, we also calculated the bias for FAPAR green , FAPAR woody , and FAPAR noWAI using the LESS-simulated FAPAR values.Table 2 summarizes the retrieval accuracy for these three methods.The results show that the TriLay method gave the best FAPAR retrieval results, having the smallest bias and RMSE values (RMSE = 0.064 and 0.025, and bias = −6.02%,−4.04% for FAPAR green under black-sky and white-sky conditions, respectively).

Temporal Variations in Different FAPAR Products
Using the TriLay approach, black-sky and white-sky products, including FAPAR green , FAPAR woody , and FAPAR canopy , and also FAPAR without woody components (FAPAR noWAI ), were generated for tile h10v05 in 2017.The black-sky FAPAR products were determined according to the SZA values at 10:30 am local time (the overpass time of the Terra satellite).
To investigate the seasonal variations in FAPAR green and the other FAPAR products, only the mean black-sky FAPAR values were calculated for the different forest vegetation types, as shown in Figure 7.The proportion of black-sky FAPAR green in FAPAR canopy and FAPAR noWAI , and also the bias between the black-sky FAPAR green , FAPAR canopy , and FAPAR noWAI were calculated for different periods during 2017 in order to analyze the differences between different black-sky FAPAR products, as well as to quantify the contribution of the woody components for several forest vegetation types (deciduous broadleaf forest, deciduous needleleaf forest, evergreen broadleaf forest, and evergreen needleleaf forest, referred to as DBF, DNF, EBF, and ENF, respectively).The results are illustrated in Figure 8 and Table 3.
In general, the black-sky FAPAR green and FAPAR noWAI exhibit typical seasonal variations for the selected forest types, with low values during the early and late period (January-February-March (JFM) and October-November-December (OND)) and high values during the peak growth stage (July-August-September (JAS)).The black-sky FAPAR canopy is obviously higher than the black-sky FAPAR green during the whole year and has a much smaller seasonal variation.The black-sky FAPAR woody behaves in the opposite way; for deciduous seasons (JFM and OND), FAPAR woody is close to or higher than its value during the peak growth stage (JAS) because the black-sky FAPAR increases with increasing SZA value.(The mean SZA at 10:30 am varies from 22.98 • (22 December) to 62.25 • (22 June) within tile h10v05.) The black-sky FAPAR green is about 52.59% and 60.60% of the black-sky FAPAR canopy for deciduous broadleaf forest during JFM and OND, respectively; the corresponding figures for deciduous needleleaf forest are only about 23.86% and 35.75%.During the peak growth stage (JAS), the black-sky FAPAR green is about 93.36%, 75.02%, 90.93% and 87.14% of FAPAR canopy for DBF, DNF, EBF, and ENF, respectively.
There are also small discrepancies between FAPAR noWAI and FAPAR green (Figure 8).In particular, for deciduous needleleaf forests, the black-sky FAPAR noWAI is overestimated by 38.30% and 28.46% during the early and late stages of the year (JFM and OND).For evergreen forests, the difference can be neglected as there is only a very slight underestimation of 0.68% to 2.39% during the whole year.

Uncertainty in Determining WAI
Currently, methods for measuring the woody area index (WAI) and woody-to-total area ratio can be classified into direct methods (e.g., destructive sampling) and indirect methods [66].However, both the direct and indirect methods can only be applied at small scales and to a limited range of vegetation types.In this study, we used a constant value of the woody-to-total area ratio for each forest type and derived the WAI from the PAI value for the peak growth stage.We obtained woody-to-total area ratios for different forest types (ENF, EBF, DBF, and DNF) from an extensive literature review [58,[64][65][66][67]-the statistical metrics, including the mean, standard deviation, and coefficient of variation, are shown in Table 4.These results show that the woody-to-total area ratios for different forest types vary from 0.158 to 0.3.The variation in this ratio within each of the forest types is also large-10.00%to 82.22%.Therefore, the woody-to-total area ratio not only varies with the forest type but also changes a lot for each individual forest type.This means that there is definitely some uncertainty due to these factors.Even so, Zou et al. [66] showed that the woody-to-total area ratio is relatively stable for the same forest stand and thus the assumption of a fixed woody-to-total area ratio for each forest type is reasonable.Differences in in situ measurement methods can also contribute to the variation in the ratio for a given forest type.The retrieval of FAPAR green can thus be improved if an accurate woody-to-total area ratio dataset is available.The canopy directional transmittance can be determined using the gap fraction model [16,53,54].Splitting the canopy into green components and woody components requires that the transmittances are also calculated separately.For the green components, the leaf absorptance at the VIS band varies from 0.79 to 0.94 as the chlorophyll content varies (from 20 to 100 µg/cm 2 under natural conditions), according to figures obtained using the PROSPECT-5 model [56,57].Also, because no remote sensing leaf chlorophyll content product is available, the extinction coefficient for leaves (k 1 ) was assumed to have a fixed value of 0.88 (corresponding to a chlorophyll content of 35 µg/cm 2 [15]).For the extinction coefficient of woody components (k 2 ), we also used a fixed value of 0.91, based on the data simulated by LESS (a needleleaf forest scene).However, Suwa et al. [68] reported a smaller extinction coefficient of 0.77 for brighter woody stems.Therefore, the extinction coefficient of woody components may vary with forest type, and thus the use of a fixed value for k 2 is also a source of error.At present, it is still very challenging to determine extinction coefficients for different forest canopy types.
In addition, a fixed value of Albedo pure was used in the TriLay model for woody vegetation types (i.e., ENF, EBF, DBF, and DNF), and was approximated based on the dense vegetation [15].According to the statistical results obtained from the MCD43A3 albedo product by Liu et al. [15], the visible albedo of dense vegetation with a "saturated" LAI value (e.g., LAI = 6) is very low and stable, with a mean value of 0.025 and a variance of 0.007 for white-sky condition, and a mean value of 0.020 and a variance of 0.006 for black-sky conditions.Therefore, the use of the prior VIS albedo values for "pure" vegetation may introduce a very small error, but it should be negligible [14].

Setting the Clumping Index for Photosynthetic and Woody Components
The clumping index characterizes the grouping of foliage within distinct canopy structures (such as tree crowns, shrubs, and row crops) relative to a random spatial distribution of leaves and is an important structural parameter for plant canopies that can influence canopy radiation regimes [49].In our TriLay approach, the clumping index for woody components was assumed to be the same as for green leaves, and so is another definite source of error [69].However, it is currently challenging to obtain the clumping effects of woody components within forest canopies [62].Furthermore, Chen et al. [61] indicated that the clumping of shoots in branches has a similar effect to the clumping of leaves within shoots.What's more, Zou et al. [62] found that the differences between the estimated CI for canopy and woody components was below 6% at the zenithal ranges of 0 • -75 • , and the difference was only 2% in the range of 30 • -60 • , which is quite small at most medium zenithal ranges thus can represent most actual conditions.Based on these results and the unavailability of CI datasets for forest canopies at large scales, we directly used the value of CI from Jiao et al. [49] to describe the clumping effect for both leaves and woody components.Therefore, the use of the same clumping index for both leaves and woody components is reasonable and credibe enough.

Conclusions
In this paper, a triple-source leaf-wood-soil layer (TriLay) method for separating FAPAR green and FAPAR woody from FAPAR canopy using the MODIS LAI, land cover, and non-linear spectral mixture model (NSM)-retrieved soil albedo [15] together with global CI products [49] was proposed.
According to the validation carried out using LESS-simulated FAPAR values, the TriLay FAPAR green was more accurate (R 2 = 0.937, RMSE = 0.064 and bias = −6.02%for black-sky conditions; R 2 = 0.997, RMSE = 0.025 and bias = −4.04% for white-sky conditions) than the traditional linear method (R 2 = 0.979, RMSE = 0.114 and bias = −18.04%for black-sky conditions; R 2 = 0.996, RMSE = 0.106 and bias = −16.93%for white-sky conditions), and also more accurate than FAPAR obtained without the consideration of woody components (FAPAR noWAI ) (R 2 = 0.920, RMSE = 0.071 and bias = −7.14%for black-sky conditions; R 2 = 0.999, RMSE = 0.043 and bias = −6.41%for white-sky conditions).A comparison of the results for black-sky FAPAR green against FAPAR noWAI and FAPAR canopy showed that the discrepancies between the black-sky FAPAR green and other FAPAR products could not be ignored for forest types.In particular, for deciduous needleleaf forest, the black-sky FAPAR green contributed only about 23.86% and 35.75% of FAPAR canopy during the early and late stages (JFM and OND) of the year, respectively, and 75.02% during the peak growth stage (JAS).There were also smaller discrepancies between the black-sky FAPAR noWAI and FAPAR green .For deciduous needleleaf forests, in particular, the black-sky FAPAR noWAI was overestimated by 38.30% and 28.46%, respectively, during the early and late stages of the year (JFM and OND).
Overall, this study provides a new method for partitioning FAPAR canopy into FAPAR green and FAPAR woody for forest types and will improve the understanding of energy exchange within the canopy.In addition, the exclusion of the contribution of woody components may certainly improve the accuracy of the FAPAR green estimates for forest types, which is significant in terms of the better modeling of vegetation photosynthesis.

Figure 2 .
Figure 2. Flowchart illustrating the Triple-source leaf-wood-soil layer (TriLay) method for estimating the fraction of radiation absorbed by photosynthetic components (FAPAR green ) and the fraction of radiation absorbed by woody components (FAPAR woody ).

Figure 3 .
Figure 3. Illustration of the triple-source leaf-wood-soil layer model.

Figure 5 .
Figure 5. Validation of FAPAR green and FAPAR woody estimated by the linear method using the LESS-simulated values of FAPAR: (a-b) black-sky FAPAR green and FAPAR woody against the LESS-simulated FAPAR green ; (c-d) white-sky FAPAR green and FAPAR woody against the LESS-simulated FAPAR woody .

Figure 8 .
Figure 8. Bias between the black-sky FAPAR green and other black-sky FAPAR products within tile h10v05 during different periods in 2017: (a) bias between FAPAR green and FAPAR noWAI ; (b) bias between FAPAR green and FAPAR canopy .JFM, AMJ, JAS, and OND represent the four seasons January to March, April to June, July to September and October to December, respectively.

Table 1 .
The main input parameters used in the LargE-Scale remote sensing data and image simulation (LESS) model simulations.

Table 2 .
Retrieval accuracy of FAPAR green , FAPAR woody , and FAPAR noWAI validated using the LESS simulations.

Table 3 .
The ratios of black-sky FAPAR green to black-sky FAPAR canopy (R canopy ) and to FAPAR noWAI (R noWAI ) for selected forest types for different periods of 2017.

Table 4 .
Statistical details of prior woody-to-total area ratios.Uncertainty Caused by the Use of Fixed Values of the Extinction Coefficients and Albedo pure