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

Abstract

The fraction of absorbed photosynthetically active radiation (FAPAR) is generally divided into the fraction of radiation absorbed by the photosynthetic components ($FAPA{R}_{green}$) and the fraction of radiation absorbed by the non-photosynthetic components ($FAPA{R}_{woody}$) 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 $FAPA{R}_{canopy}$ into $FAPA{R}_{green}$ and $FAPA{R}_{woody}$ 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 $FAPA{R}_{green}$ had higher retrieval accuracy, as well as a significantly lower bias (R² = 0.937, Root Mean Square Error (RMSE) = 0.064, and bias = −6.02% for black-sky conditions; R² = 0.997, RMSE = 0.025 and bias = −4.04% for white-sky conditions) compared to the traditional linear method (R² = 0.979, RMSE = 0.114, and bias = −18.04% for black-sky conditions; R² = 0.996, RMSE = 0.106 and bias = −16.93% for white-sky conditions). For FAPAR that did not take account of woody components ($FAPA{R}_{noWAI}$), the corresponding results were R² = 0.920, RMSE = 0.071, and bias = −7.14% for black-sky conditions, and R² = 0.999, RMSE = 0.043, and bias = −6.41% for white-sky conditions. Finally, the dynamic $FAPA{R}_{green}$, $FAPA{R}_{woody}$, $FAPA{R}_{canopy}$ and $FAPA{R}_{noWAI}$ 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 $FAPA{R}_{green}$ against those for $FAPA{R}_{noWAI}$ and $FAPA{R}_{canopy}$ showed that the discrepancy between $FAPA{R}_{green}$ and other FAPAR products for forest vegetation types could not be ignored. For deciduous needleleaf forest, in particular, the black-sky $FAPA{R}_{green}$ was found to contribute only about 23.86% and 35.75% of $FAPA{R}_{canopy}$ 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 $FAPA{R}_{noWAI}$ 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 $FAPA{R}_{green}$ from $FAPA{R}_{canopy}$, which is of great importance for a better understanding of the energy exchange within the canopy.

1. 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].

Recently, several global FAPAR products have become available, including the Moderate Resolution Imaging Spectroradiometer (MODIS) [20,21], Energy Balance Residual (EBR) [15], Multi-angle Imaging SpectroRadiometer (MISR) [11], CYCLOPES [22], GLOBCARBON [23], Global Land Surface Satellite (GLASS) [14], the Medium Resolution Imaging Spectrometer (MERIS) [24], Joint Research Center Two Stream Inversion Package (JRC-TIP) [25], and European Space Agency (ESA) products [26]. These global FAPAR products have been widely validated, with reported errors varying from 0.08 to 0.23 [14,27,28,29,30,31,32]. However, most of the global FAPAR products do not consider the effect of non-photosynthetic components in the radiative transfer process, which introduces errors, especially for forest types. Moreover, many researchers have used the fraction of radiation absorbed by photosynthetic components ($FAPA{R}_{green}$) instead of the fraction of radiation absorbed by the canopy ($FAPA{R}_{canopy}$) to monitor and estimate the light use efficiency (LUE), radiation use efficiency (RUE), and productivity at different temporal scales [33,34,35,36].

The forest vegetation ecosystem plays an important role in the global ecosystem. However, quantifying the temporal variation in $FAPA{R}_{green}$ for a forest ecosystem represents an important challenge for remote sensing and ecology researchers as it is extremely difficult to measure $FAPA{R}_{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 $FAPA{R}_{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 $FAPA{R}_{green}$ from remote sensing data. Hall et al. [38] estimated $FAPA{R}_{green}$ using a simple linear relationship between $FAPA{R}_{canopy}$ and $LA{I}_{green}/LA{I}_{total}$ ($LA{I}_{total}$ denotes the total leaf area index including green and senescent leaves, while $LA{I}_{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 $FAPA{R}_{green}$ and $FAPA{R}_{canopy}$ using the forward simulation approach. However, $FAPA{R}_{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 $FAPA{R}_{canopy}$ into photosynthetically active green components ($FAPA{R}_{green}$) and non-photosynthetic active components using the ratio $LA{I}_{green}/LA{I}_{total}$ for maize and soybeans. The relationship between vegetation indices and $FAPA{R}_{green}$ was also used to retrieve $FAPA{R}_{green}$ [41,42]. Nevertheless, to date, the current $FAPA{R}_{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 $FAPA{R}_{canopy}$ into $FAPA{R}_{green}$ and $FAPA{R}_{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. $FAPA{R}_{canopy}$ is first separated into the fraction of PAR absorbed by the canopy for downwelling radiation ($FAPA{R}_{canopy\downarrow }$) and the fraction of PAR absorbed by the canopy for the upwelling radiation reflected by the soil background ($FAPA{R}_{canopy\uparrow }$). Then, $FAPA{R}_{canopy\downarrow }$ and $FAPA{R}_{canopy\uparrow }$ are further split into the fraction of radiation absorbed by photosynthetic components ($FAPA{R}_{green}$) and that absorbed by non-photosynthetic components ($FAPA{R}_{woody}$) using the TriLay model. Finally, the $FAPA{R}_{green}$, $FAPA{R}_{woody}$, and $FAPA{R}_{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 $FAPA{R}_{green}$ and $FAPA{R}_{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 $FAPA{R}_{green}$ estimates, which is important for better modeling of vegetation photosynthesis.

2. Materials and Methods

2.1. Satellite Datasets

In order to produce $FAPA{R}_{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 $FAPA{R}_{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 $FAPA{R}_{green}$ and other FAPAR products because there were abundant forest vegetation types within this MODIS tile.

2.1.1. 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].

2.1.2. 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.

2.1.3. 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].

2.1.4. 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.

2.2. Data Simulated by the LESS Model

To quantitatively evaluate the performance of the proposed TriLay approach for estimating $FAPA{R}_{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²) 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 $FAPA{R}_{canopy}$, $FAPA{R}_{green}$ and $FAPA{R}_{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. The main input parameters used in the LargE-Scale remote sensing data and image simulation (LESS) model simulations.

Parameter	Definition	Units	Range or Values
Canopy
LAI	leaf area index	m2/m2	1.31–8.69
WAI	woody area index	m2/m2	1.65
Leaf layer
Reflectance		—	0.041–0.205
Transmittance		—	0.001–0.286
Soil layer
Reflectance		—	0.001–0.134
Woody layer
Reflectance		—	0.069–0.237
Imaging Geometry
SZA	sun zenith angle	degrees	0, 10, 20, 30, 40, 50, 60, 70, 80
$Rati{o}_{Sky}$	ratio of diffuse light	—	0, 1

. Finally, a total of 72 black-sky simulations and eight white-sky simulations were achieved for the different conditions giving.

2.3. Algorithms for Estimating Global $FAPA{R}_{green}$ and $FAPA{R}_{woody}$ Datasets

To split the fraction of radiation absorbed by photosynthetic components ($FAPA{R}_{green}$) from $FAPA{R}_{canopy}$, a novel triple-source leaf–wood–soil layer model was proposed for generating global $FAPA{R}_{green}$ products for forest vegetation types. Figure 2 is a flowchart of the process used to retrieve global $FAPA{R}_{green}$ products.

First, $FAPA{R}_{canopy}$ is split into two parts: the fraction of PAR absorbed by the canopy for the downwelling radiation ($FAPA{R}_{canopy\downarrow }$) and that absorbed by the canopy for the upwelling radiation reflected by the soil background ($FAPA{R}_{canopy\uparrow }$). Then, $FAPA{R}_{canopy\downarrow }$ and $FAPA{R}_{canopy\uparrow }$ are further split into the fraction of radiation absorbed by photosynthetic components ($FAPA{R}_{green}$) and that absorbed by non-photosynthetic components (e.g., branches and stems, hereafter called $FAPA{R}_{woody}$). Finally, $FAPA{R}_{green}$ and $FAPA{R}_{woody}$ can be calculated separately using the TriLay approach.

2.3.1. 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, $FAPA{R}_{canopy}$, can be separated into the fraction of PAR absorbed by the canopy for the downwelling radiation ($FAPA{R}_{canopy\downarrow }$) and that for the upwelling radiation reflected from the soil background ($FAPA{R}_{canopy\uparrow }$):

$$\begin{array}{c}FAPA{R}_{canopy}=FAPA{R}_{canopy\downarrow }+FAPA{R}_{canopy\uparrow }\end{array}$$

(1)

where $FAPA{R}_{canopy\downarrow }$ describes the fraction of downwelling radiation absorbed within the canopy assuming the soil background is dark, and $FAPA{R}_{canopy\uparrow }$ 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 $FAPA{R}_{canopy\downarrow }$ is equal to the $FAPA{R}_{canopy}$.

Therefore, $FAPA{R}_{canopy\downarrow }$ can be given by:

$$\begin{array}{c}FAPA{R}_{canopy}=\left(1-{\tau }_{PAI}\right)\ast \left(1-Albed{o}_{pure}\ast FVC\right)\end{array}$$

(2)

where ${\tau }_{PAI}$ is the transmittance of the whole canopy, which contains green and woody parts, $FVC$ is the fraction of vegetation cover, and $Albed{o}_{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], $Albed{o}_{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).

${\tau }_{PAI}$ can be calculated as the product of ${\tau }_{LAI}$ and ${\tau }_{WAI}$, and the directional transmittance can be determined using the gap fraction model [16,53,54]:

$$\begin{array}{c}{\tau }_{PAI}={\tau }_{LAI}\times {\tau }_{WAI}\end{array}$$

(3)

$$\begin{array}{c}{\tau }_{LAI}=e^{-k_1\times G\left(\theta \right)\times CI\times LAI/\cos \left(\theta \right)}\end{array}$$

(4)

$$\begin{array}{c}{\tau }_{WAI}=e^{-k_2\times G\left(\theta \right)\times CI\times WAI/\cos \left(\theta \right)}\end{array}$$

(5)

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 $\theta$ 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\left(\theta \right)$ is the projection of the unit foliage area on the plane perpendicular to the solar incident direction, $\theta$. For green leaves, $G\left(\theta \right)$ 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:

$$\begin{array}{c}FVC=1-e^{-G\left(\theta \right)\times LAI\times CI}\end{array}$$

(6)

Similarly,$FAPA{R}_{canopy\uparrow }$ can be calculated as follows:

$$\begin{array}{c}FAPA{R}_{canopy\uparrow }=\left(1-{\tau }_{PAI}\right)\ast \left(1-Albed{o}_{pure}\ast \mathrm{FVC}\right)\ast {\tau }_{PAI}^{ws}\ast Albed{o}_{soil}\end{array}$$

(7)

where $Albed{o}_{soil}$ is the soil visible albedo, which can be generated by the NSM model proposed by Liu et al. [15]. ${\tau }_{PAI}^{ws}$ is the canopy transmittance under white-sky conditions, which can be calculated as the product of the white-sky transmittance of leaves (${\tau }_{LAI}^{ws}$) and the white-sky transmittance of woody components (${\tau }_{WAI}^{ws}$):

$$\begin{array}{c}{\tau }_{LAI}^{ws}=2\times {{\displaystyle \int }}_0^{\frac{\pi }{2}}\left(e^{-k_1\times G\left(\theta \right)\times LAI\times CI/\cos \left(\theta \right)}\right)\times \sin \left(\theta \right)\times \cos \left(\theta \right)d\theta \end{array}$$

(8)

$$\begin{array}{c}{\tau }_{WAI}^{ws}=2\times {{\displaystyle \int }}_0^{\frac{\pi }{2}}\left(e^{-k_2\times G\left(\theta \right)\times WAI\times CI/\cos \left(\theta \right)}\right)\times \sin \left(\theta \right)\times \cos \left(\theta \right)d\theta \end{array}$$

(9)

$$\begin{array}{c}{\tau }_{PAI}^{ws}={\tau }_{LAI}^{ws}\times {\tau }_{WAI}^{ws}\end{array}$$

(10)

2.3.2. 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 $FAPA{R}_{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-to- total 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:

$$\begin{array}{c}WAI=PAI\times rati{o}_{woody}=\frac{LA{I}_{max}}{1-rati{o}_{woody}}\times rati{o}_{woody}\end{array}$$

(11)

where $rati{o}_{woody}$ is the mean value of the woody-to-total area ratio for various forest types, as given in the literature [58,64,65]. $LA{I}_{max}$ is the maximum value of LAI within a given year; this was acquired from MCD15A2H LAI products.

2.3.3. Separating $FAPA{R}_{green}$ and $FAPA{R}_{woody}$ from $FAPA{R}_{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 $FAPA{R}_{canopy\downarrow }$ is split into $FAPA{R}_{green\downarrow }$ and $FAPA{R}_{woody\downarrow }$:

$$\begin{array}{c}FAPA{R}_{canopy\downarrow }=FAPA{R}_{green\downarrow }+FAPA{R}_{woody\downarrow }\end{array}$$

(12)

where $FAPA{R}_{green\downarrow }$ and $FAPA{R}_{woody\downarrow }$ can be obtained as

$$\begin{array}{c}FAPA{R}_{green\downarrow }=FAPA{R}_{canopy\downarrow }\times rati{o}_{green}\times w_{1\downarrow }\end{array}$$

(13)

$$\begin{array}{c}FAPA{R}_{woody\downarrow }=FAPA{R}_{canopy\downarrow }\times rati{o}_{woody}\times w_{2\downarrow }\end{array}$$

(14)

where $rati{o}_{green}$ and $rati{o}_{woody}$ are the ratio of leaf area index to plant area index and woody area index to plant area index, respectively. $w_{1\downarrow }$ and $w_{2\downarrow }$ 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:

$$\begin{array}{c} rati{o}_{green}=\frac{LAI}{PAI}\end{array}$$

(15)

$$\begin{array}{c}rati{o}_{woody}=\frac{WAI}{PAI}\end{array}$$

(16)

$$\begin{array}{c}{w}_{2\downarrow }=w_{1\downarrow }\times {\tau }_{LAI}\end{array}$$

(17)

$$\begin{array}{c}{w}_{1\downarrow }\times rati{o}_{green}+w_{2\downarrow }\times rati{o}_{woody}=1\end{array}$$

(18)

$FAPA{R}_{green\downarrow }$ and $FAPA{R}_{woody\downarrow }$ can then be obtained by solving equations (13)–(18):

$$\begin{array}{c}FAPA{R}_{green\downarrow }=\frac{rati{o}_{green}\times FAPA{R}_{canopy\downarrow }}{rati{o}_{green}+{\tau }_{LAI}\times rati{o}_{woody}}\end{array}$$

(19)

$$\begin{array}{c}FAPA{R}_{woody\downarrow }=\frac{rati{o}_{woody}\times FAPA{R}_{canopy\downarrow }\times {\tau }_{LAI}}{rati{o}_{green}+{\tau }_{LAI}\times rati{o}_{woody}}\end{array}$$

(20)

Similarly, $FAPA{R}_{green\uparrow }$ and $FAPA{R}_{woody\uparrow }$ can also be acquired:

$$\begin{array}{c}FAPA{R}_{green\uparrow }=\frac{rati{o}_{green}\times FAPA{R}_{canopy\uparrow }\times {\tau }_{WAI}}{rati{o}_{woody}+{\tau }_{WAI}\times rati{o}_{green}}\end{array}$$

(21)

$$\begin{array}{c}FAPA{R}_{woody\uparrow }=\frac{rati{o}_{woody}\times FAPA{R}_{canopy\uparrow }}{rati{o}_{woody}+{\tau }_{WAI}\times rati{o}_{green}}\end{array}$$

(22)

Finally, $FAPA{R}_{green}$ and $FAPA{R}_{woody}$ can be calculated as:

$$\begin{array}{c}FAPA{R}_{green}=FAPA{R}_{green\uparrow }+FAPA{R}_{green\downarrow }\end{array}$$

(23)

$$\begin{array}{c}FAPA{R}_{woody}=FAPA{R}_{woody\uparrow }+FAPA{R}_{woody\downarrow }\end{array}$$

(24)

3. Results

3.1. Validation of the TriLay Method using Simulations made by the LESS Model

It is too challenging to obtain in-situ measurements of $FAPA{R}_{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 $FAPA{R}_{canopy}$, $FAPA{R}_{green}$, and $FAPA{R}_{woody}$ together with eight white-sky simulations were available for validation.

Figure 4 illustrates the validation results for $FAPA{R}_{canopy}$, $FAPA{R}_{green}$, and $FAPA{R}_{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 $FAPA{R}_{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² values are 0.945 and 0.999. For $FAPA{R}_{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 $FAPA{R}_{woody}$ can also be accurately estimated: the RMSE and R² 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 $FAPA{R}_{canopy}$, $FAPA{R}_{green}$, and $FAPA{R}_{woody}$ for forest land cover types.

Furthermore, as illustrated in Figure 4b, there is a slight underestimation for $FAPA{R}_{green}$ at smaller SZAs (0°–60°) and a slight overestimation for larger SZAs (70°–80°). From Figure 4c, it can be seen that $FAPA{R}_{woody}$ is also slightly underestimated at smaller SZAs (0°–40°) and overestimated at larger angles (50°–80°).

3.2. Comparison of Different Methods using the LESS Simulations

Hall et al. [38] estimated $FAPA{R}_{green}$ based on a simple linear relationship between $FAPA{R}_{canopy}$ and $FAPA{R}_{green}$. They used the ratio $LA{I}_{green}/LA{I}_{total}$ to determine $FAPA{R}_{green}$:

$$\begin{array}{c}FAPA{R}_{green}=FAPA{R}_{canopy}\times \frac{LA{I}_{green}}{totalLAI}\end{array}$$

(25)

where total LAI is the PAI mentioned above. In order to test the accuracy of the linear mixture method, we also validated the $FAPA{R}_{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 $FAPA{R}_{green}$ and a huge overestimation for $FAPA{R}_{woody}$. Although the FAPARs retrieved using the linear mixture method are highly correlated with the LESS simulation, with R² values of 0.979, 0.996 for $FAPA{R}_{green}$, and 0.934, 0.985 for $FAPA{R}_{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 $FAPA{R}_{green}$ as against 0.064, 0.025, and 0.113, 0.106 for $FAPA{R}_{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 ($FAPA{R}_{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² 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² 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 $FAPA{R}_{green}$, $FAPA{R}_{woody}$, and $FAPA{R}_{noWAI}$ using the LESS-simulated FAPAR values.

Table 2. Retrieval accuracy of $FAPA{R}_{green}$, $FAPA{R}_{woody}$, and $FAPA{R}_{noWAI}$ validated using the LESS simulations.

(a) For $FAPA{R}_{green}$ Products
$\mathit{F}\mathit{A}\mathit{P}\mathit{A}{\mathit{R}}_{\mathit{g}\mathit{r}\mathit{e}\mathit{e}\mathit{n}}$	TriLay		Linear		noWAI
	Black-Sky	White-Sky	Black-Sky	White-Sky	Black-Sky	White-Sky
R2	0.937	0.997	0.979	0.996	0.920	0.999
RMSE	0.064	0.025	0.114	0.106	0.071	0.043
Bias	−6.02%	−4.04%	−18.04%	−16.93%	−7.14%	−6.41%
(b) For $FAPA{R}_{woody}$ products
$\mathit{F}\mathit{A}\mathit{P}\mathit{A}{\mathit{R}}_{\mathit{w}\mathit{o}\mathit{o}\mathit{d}\mathit{y}}$	TriLay		Linear
	Black-Sky	White-Sky	Black-Sky	White-Sky
R2	0.709	0.992	0.934	0.985
RMSE	0.042	0.014	0.113	0.106
Bias	6.87%	−4.64%	153.84%	123.47%

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 $FAPA{R}_{green}$ under black-sky and white-sky conditions, respectively).

3.3. Temporal Variations in Different FAPAR Products

Using the TriLay approach, black-sky and white-sky products, including $FAPA{R}_{green}$, $FAPA{R}_{woody}$, and $FAPA{R}_{canopy}$, and also FAPAR without woody components ($FAPA{R}_{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 $FAPA{R}_{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 $FAPA{R}_{green}$ in $FAPA{R}_{canopy}$ and $FAPA{R}_{noWAI}$, and also the bias between the black-sky $FAPA{R}_{green}$, $FAPA{R}_{canopy},$ and $FAPA{R}_{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. The ratios of black-sky $FAPA{R}_{green}$ to black-sky $FAPA{R}_{canopy}$ ($R_{canopy}$) and to $FAPA{R}_{noWAI}$ ($R_{noWAI}$) for selected forest types for different periods of 2017.

Period of Year	DBF		DNF		EBF		ENF
	${\mathit{R}}_{\mathit{c}\mathit{a}\mathit{n}\mathit{o}\mathit{p}\mathit{y}}$	${\mathit{R}}_{\mathit{n}\mathit{o}\mathit{W}\mathit{A}\mathit{I}}$	${\mathit{R}}_{\mathit{c}\mathit{a}\mathit{n}\mathit{o}\mathit{p}\mathit{y}}$	${\mathit{R}}_{\mathit{n}\mathit{o}\mathit{W}\mathit{A}\mathit{I}}$	${\mathit{R}}_{\mathit{c}\mathit{a}\mathit{n}\mathit{o}\mathit{p}\mathit{y}}$	${\mathit{R}}_{\mathit{n}\mathit{o}\mathit{W}\mathit{A}\mathit{I}}$	${\mathit{R}}_{\mathit{c}\mathit{a}\mathit{n}\mathit{o}\mathit{p}\mathit{y}}$	${\mathit{R}}_{\mathit{n}\mathit{o}\mathit{W}\mathit{A}\mathit{I}}$
	(%)		(%)		(%)		(%)
JFM	52.59	93.14	23.86	73.90	82.94	101.01	74.55	100.72
AMJ	90.74	101.93	64.13	96.32	91.03	102.27	86.85	102.46
JAS	93.36	102.24	75.02	99.19	90.93	102.19	87.14	102.46
OND	60.60	95.50	35.75	81.65	84.54	101.20	78.00	101.19

JFM, AMJ, JAS, and OND represent the four seasons January to March, April to June, July to September and October to December, respectively.

.

In general, the black-sky $FAPA{R}_{green}$ and $FAPA{R}_{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 $FAPA{R}_{canopy}$ is obviously higher than the black-sky $FAPA{R}_{green}$ during the whole year and has a much smaller seasonal variation. The black-sky $FAPA{R}_{woody}$ behaves in the opposite way; for deciduous seasons (JFM and OND), $FAPA{R}_{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 $FAPA{R}_{green}$ is about 52.59% and 60.60% of the black-sky $FAPA{R}_{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 $FAPA{R}_{green}$ is about 93.36%, 75.02%, 90.93% and 87.14% of $FAPA{R}_{canopy}$ for DBF, DNF, EBF, and ENF, respectively.

There are also small discrepancies between $FAPA{R}_{noWAI}$ and $FAPA{R}_{green}$ (Figure 8). In particular, for deciduous needleleaf forests, the black-sky $FAPA{R}_{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.

4. Discussion

4.1. 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. Statistical details of prior woody-to-total area ratios.

Forest Vegetation Type	ENF	EBF	DNF	DBF
number of samples	35	8	3	4
mean value	0.185	0.18	0.3	0.158
standard deviation	0.062	0.148	0.03	0.101
coefficient of variation	33.51%	82.22%	10.00%	63.92%

. 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 $FAPA{R}_{green}$ can thus be improved if an accurate woody-to-total area ratio dataset is available.

4.2. Uncertainty Caused by the Use of Fixed Values of the Extinction Coefficients and $Albed{o}_{pure}$

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² 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² [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 $Albed{o}_{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].

4.3. 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.

5. Conclusions

In this paper, a triple-source leaf–wood–soil layer (TriLay) method for separating $FAPA{R}_{green}$ and $FAPA{R}_{woody}$ from $FAPA{R}_{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 $FAPA{R}_{green}$ was more accurate (R² = 0.937, RMSE = 0.064 and bias = −6.02% for black-sky conditions; R² = 0.997, RMSE = 0.025 and bias = −4.04% for white-sky conditions) than the traditional linear method (R² = 0.979, RMSE = 0.114 and bias = −18.04% for black-sky conditions; R² = 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 ($FAPA{R}_{noWAI}$) (R² = 0.920, RMSE = 0.071 and bias = −7.14% for black-sky conditions; R² = 0.999, RMSE = 0.043 and bias = −6.41% for white-sky conditions). A comparison of the results for black-sky $FAPA{R}_{green}$ against $FAPA{R}_{noWAI}$ and $FAPA{R}_{canopy}$ showed that the discrepancies between the black-sky $FAPA{R}_{green}$ and other FAPAR products could not be ignored for forest types. In particular, for deciduous needleleaf forest, the black-sky $FAPA{R}_{green}$ contributed only about 23.86% and 35.75% of $FAPA{R}_{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 $FAPA{R}_{noWAI}$ and $FAPA{R}_{green}$. For deciduous needleleaf forests, in particular, the black-sky $FAPA{R}_{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 $FAPA{R}_{canopy}$ into $FAPA{R}_{green}$ and $FAPA{R}_{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 $FAPA{R}_{green}$ estimates for forest types, which is significant in terms of the better modeling of vegetation photosynthesis.