Potential Investigation of Linking PROSAIL with the Ross-Li BRDF Model for Vegetation Characterization

Abstract

Methods that link different models for investigating the retrieval of canopy biophysical/structural variables have been substantially adopted in the remote sensing community. To retrieve global biophysical parameters from multiangle data, the kernel-driven bidirectional reflectance distribution function (BRDF) model has been widely applied to satellite multiangle observations to model (interpolate/extrapolate) the bidirectional reflectance factor (BRF) in an arbitrary direction of viewing and solar geometries. Such modeled BRFs, as an essential information source, are then input into an inversion procedure that is devised through a large number of simulation analyses from some widely used physical models that can generalize such an inversion relationship between the BRFs (or their simple algebraic composite) and the biophysical/structural parameter. Therefore, evaluation of such a link between physical models and kernel-driven models contributes to the development of such inversion procedures to accurately retrieve vegetation properties, particularly based on the operational global BRDF parameters derived from satellite multiangle observations (e.g., MODIS). In this study, the main objective is to investigate the potential for linking a popular physical model (PROSAIL) with the widely used kernel-driven Ross-Li models. To do this, the BRFs and albedo are generated by the physical PROSAIL in a forward model, and then the simulated BRFs are input into the kernel-driven BRDF model for retrieval of the BRFs and albedo in the same viewing and solar geometries. To further strengthen such an investigation, a variety of field-measured multiangle reflectances have also been used to investigate the potential for linking these two models. For simulated BRFs generated by the PROSAIL model at 659 and 865 nm, the two models are generally comparable to each other, and the resultant root mean square errors (RMSEs) are 0.0092 and 0.0355, respectively, although some discrepancy in the simulated BRFs can be found at large average leaf angle (ALA) values. Unsurprisingly, albedos generated by the method are quite consistent, and 99.98% and 97.99% of the simulated white sky albedo (WSA) has a divergence less than 0.02. For the field measurements, the kernel-driven model presents somewhat better model-observation congruence than the PROSAIL model. The results show that these models have an overall good consistency for both field-measured and model-simulated BRFs. Therefore, there is potential for linking these two models for looking into the retrieval of canopy biophysical/structural variables through a simulation method, particularly from the current archive of the global routine MODIS BRDF parameters that were produced by the kernel-driven BRDF model; however, erectophile vegetation must be further examined.

1. Introduction

Linking different models for various applications in remote sensing is an important method. Typical examples arise from linking leaf optical properties models with canopy bidirectional reflectance models to characterize the entire process of light absorption and scattering of the vegetation canopy. The benefit of this technique also makes it possible to retrieve both leaf biochemistry characteristics and canopy structures from satellite data [1]. By linking the optical Properties Spectra (PROSPECT) broadleaf model and the homogeneous Scattering by Arbitrarily Inclined Leaves (SAIL) canopy reflectance (CR) model, the PROPECT+SAIL (named PROSAIL) model has been widely used to simulate canopy spectral and directional reflectance and to retrieve biophysical parameters [1,2,3]. The practicability of this model has been substantially validated using large numbers of field measurements. Specifically, at the leaf level, various leaf observations from the LOPEX, CALMIT, ANGERS and HAWAII campaigns launched from 1993 to 2007 have been used for validation [2], and at the canopy level, the model was comparable to many other 1-D and 3-D models in the famous controlled Radiation Transfer Model Intercomparison (RAMI) experiments from Phase 1 to Phase 3 [4,5,6]. This model was also used in conjunction with a kernel-driven model to evaluate the emissivity estimation in cases with insufficient field data [7], and to assess a novel Angular and Spectral Kernel-driven (ASK) model [8]. Similarly, the PROSPECT model has also been linked with the A two-layer Canopy Reflectance Model (ACRM) [9], the New Advanced Discrete (NADI) model [10] and the 3-D Discrete Anisotropic Radiative Transfer (DART) model (URL: http://teledetection.ipgp.jussieu.fr/prosail/) for various applications. More extensively, the other conifer leaf model named the LIBERTY model was linked with the 4-Scale CR model to generate the so-called 5-Scale model [11,12,13] for the enhancement of simulation ability of these models.

To implement quantitative remote sensing inversion from reflection signatures of land surface, the intrinsic anisotropy reflection characteristics of land surface play a key role, which can be characterized by Bidirectional Reflectance Distribution Function (BRDF) successfully [14]. As an essential source of BRDF data, available remotely sensed multiangle observations and BRDF products at the global scale (e.g., satellite multiangle data from the Moderate-resolution Imaging Spectroradiometer (MODIS) and Polarization and Directionality of the Earth’s Reflectances (POLDER)) contain important biophysical/structural information of the land surface, which is usually used as primary input parameters of physical BRDF models to drive inversion procedures [15,16,17,18,19,20]. Therefore, direct multiangle observations from satellite sensors are desperately needed to drive model inversion procedures for estimation of various biophysical/structural variables of the land surface. For example, the bi-directional reflectance product of MODIS (i.e., MOD09A1) has typically been used for LAI and albedo estimation through the ACRM [21] as well as the SAIL CR model, which were reported to yield high-accuracy retrievals based on such direct satellite multiangle observations [22,23,24].

However, in situations where satellite sensors cannot provide direct observations due to their insufficient BRDF sampling capacity (e.g., MODIS), or cloud contaminations or poor atmospheric conditions, a kernel-driven BRDF model provides an essential tool to derive the wanted BRFs in a given viewing and solar geometries by an interpolation or extrapolation method [18]. Such modeled BRFs are then usually input into the inversion procedure of physical model. A typical example is the accurate retrieval of the clumping index (CI) of vegetation foliage through use of the BRFs in the directions of Hotspot and Darkspot that are used to calculate the Normalized Difference between the Hotspot and Darkspot (NDHD) index [25,26]. This NDHD index was, in return, utilized to map the global CI through a simulation inversion procedure built by the 4-Scale physical CR model [26,27]. Likewise, by simulating customized directional reflectance as the training dataset for a physical process in regard to directional escape probabilities, canopy height at three typical forest regions was accurately estimated [28].

As a result, the method for linking a physical model with the kernel-driven BRDF model provides potentials to accurately retrieve canopy biophysical/structural variables, particularly from the global routine BRDF parameter products generated by the kernel-driven BRDF model. Despite its great potential mentioned above, it is also necessary to investigate the limitations of this method, because the linked models usually have their own specific parameterization method in modeling inhomogeneity of the land surface. As one of the popular physical BRDF models, the PROSAIL model is designed based on radiative transfer theory, where only the leaves are included with a random azimuth distribution [3,29], and one of intrinsic limitations in its capacity may arise from simulating heterogeneous canopies showing clumping at several scales [1]. For semi-empirical, kernel-driven models such as the RTLSR model, although both of the radiative transfer process and the geometric-optical mutual shadowing have been considered, the substantial simplification into several so-called kernels [30,31] may make the models inflexible to provide a fit to any type of reasonable shape, but rather are constrained by the predetermined BRDF shapes of these kernels [32]. Therefore, such limitations behind the model mechanism may yield inconsistent BRDF results once they are linked in terms of the input-output chain, which probably leads to uncertainties and errors in their applications.

In this study, the BRF and albedo simulations derived by linking the PROSAIL model with the kernel-driven Ross-Li models are compared for the first time to examine their consistency in BRFs. First, 20,000 sets of BRFs were produced through the simulation method by the PROSAIL model, which covered various optical and structural parameters of vegetation canopies and contained 397 angles for each data set. Then, these data were input into the kernel-driven models for calculating the BRFs in the same solar and viewing geometries. The albedos were then calculated from the BRFs of both models. Finally, we compared the two models in terms of their BRFs and albedos and examined how vegetation parameters affect such a model link. To further strengthen such an investigation, 21 sets of multiangle reflectances were also used in this study. Finally, we conclude this study with a summary of the findings.

2. Methods

2.1. The Kernel-Driven BRDF Model and Anisotropy Flat Index (AFX)

The kernel-driven model derived from the Roujean model contains three adjustable coefficients that are used for optimal BRDF fitting based on directional reflectance observations [30,31,33,34,35]. This model characterizes surface reflection as a simple linear combination of three types of scattering, and the general formula is given as Equation (1).

$$R({\theta }_i,{\theta }_v,\phi ,\lambda )=f_{iso}\left(\lambda \right)f_{vol}\left(\lambda \right)K_{vol}({\theta }_i,{\theta }_v,\phi )f_{geo}\left(\lambda \right)K_{geo}({\theta }_i,{\theta }_v,\phi ),$$

(1)

where $R$ is the surface directional reflectance, isotropic scattering is a constant (1.0), $K_{vol}$ is the volumetric scattering kernel caused by multiple scattering within canopies, and $K_{geo}$ represents the geometric-optical scattering kernel associated with single scattering among canopies. ${\theta }_i$, ${\theta }_v$, and $\phi$ are the incident zenith angle (generally referred to as the solar zenith angle (SZA)), reflected zenith angle (generally called the view zenith angle (VZA)) and relative azimuth angle, respectively. $f_{iso}$, $f_{vol}$ and $f_{geo}$ are the weights of the three kernels. The hotspot-revised Ross Thick-Chen (RTC) volume scattering kernel is shown in Equation (2), which is a correction of the Ross Thick model that is done by multiplying a hotspot factor with its first term [33].

$$K_{RTC}=\frac{(\pi /2-\xi )\cos \xi +\sin \xi }{\cos {\theta }_v+\cos {\theta }_i}\times (1+C_1{e}^{-\frac{\xi }{C_2}})-\frac{\pi }{4},$$

(2)

where C₁ and C₂ refer to the height and width of the hotspot, respectively, and the optimal values of the two hotspot parameters can be first determined based on the minimum root mean square error (RMSE) between hotspot observations and retrievals derived from the combined RTCLSR/RTCLTR model (given in Equation (17)) [33]. $\xi$ is the phase angle, which can be described as follows:

$$\cos \xi =\cos {\theta }_i\cos {\theta }_v+\sin {\theta }_i\sin {\theta }_v\cos \phi .$$

(3)

$K_{geo}$ was first proposed by Roujean [31], and several subsequent versions have been developed, including the sparse canopy-assuming Li Sparse-Reciprocal (LSR) kernel [17,30,34], Li Dense Reciprocal (LDR) kernel for a dense canopy [30] and sparse-dense transition Li Transit-Reciprocal (LTR) kernel [36], as shown in Equations (4)–(13) based on a literature review [37].

$$K_{LSR}=O({\theta }_i{}^{\prime },{\theta }_v{}^{\prime },t)-\sec {\theta }_i{}^{\prime }-\sec {\theta }_v{}^{\prime }+\frac{1}{2}(1+\cos {\xi }^{\prime })\sec {\theta }_v{}^{\prime }\sec {\theta }_i{}^{\prime },$$

(4)

$$K_{LDR}=\frac{(1+\cos {\xi }^{\prime })\sec {\theta }_v{}^{\prime }\sec {\theta }_i{}^{\prime }}{\sec {\theta }_i{}^{\prime }+\sec {\theta }_v{}^{\prime }-O({\theta }_i{}^{\prime },{\theta }_v{}^{\prime },t)}-2,$$

(5)

$$O({\theta }_i{}^{\prime },{\theta }_v{}^{\prime },t)=\frac{1}{\pi }(t-\sin t\cdot \cos t)(\sec {\theta }_i{}^{\prime }+\sec {\theta }_v{}^{\prime }),$$

(6)

$$\cos t=\frac{h}{b}\frac{\sqrt{D^2(\tan {\theta }_i{}^{\prime }\tan {\theta }_v{}^{\prime }\sin \phi )}}{\sec {\theta }_i{}^{\prime }+\sec {\theta }_v{}^{\prime }},$$

(7)

$$D=\sqrt{{\tan }^2{\theta }_i{}^{\prime }+{\tan }^2{\theta }_v{}^{\prime }-2\tan {\theta }_i{}^{\prime }\tan {\theta }_v{}^{\prime }\cos \phi },$$

(8)

$$\cos {\xi }^{\prime }=\cos {\theta }_i{}^{\prime }\cos {\theta }_v{}^{\prime }+\sin {\theta }_i{}^{\prime }\sin {\theta }_v{}^{\prime }\cos \phi ,$$

(9)

$${\theta }_i{}^{\prime }={\tan }^{-1}(\frac{b}{r}\tan {\theta }_i), \mathrm{and}$$

(10)

$${\theta }_v{}^{\prime }={\tan }^{-1}(\frac{b}{r}\tan {\theta }_v).$$

(11)

The dimensionless relative crown height $h/b$ and shape $b/r$ included in Equations (7), (10) and (11) are set to values of 2 and 1, respectively, for MODIS BRDF processing [17]. The LSR and LDR models are combined to develop the LTR model and perform accurate forecasting using Equation (12). The LSR model is applied at small viewing angles, and the LDR model is used at large VZAs when the critical value of B is greater than 2 [36]. Consequently, the LTR model retains the robust fitting ability of the widely used LSR model and improves the extrapolation process at large zenith angles. Therefore, to identify the model that is ideal for model coupling, the widely used LSR model, which is used as the operational algorithm of the MODIS BRDF product, and the optimal LTR model were analyzed in this study.

$$K_{LTR}=\{\begin{array}{l}{K}_{LSR},B\le 2\\ L_{LDR}=\frac{2}{B}{K}_{LSR},B>2\end{array},$$

(12)

where

$$B=\sec {\theta }_i{}^{\prime }+\sec {\theta }_v{}^{\prime }-O({\theta }_i{}^{\prime },{\theta }_v{}^{\prime },\phi ).$$

(13)

Based on Equation (1), we retrieved the three optimal kernel coefficients of multiangle observations using the least squares method, where the three coefficients were constrained to be non-negative. Then, the all-directional reflectance was calculated for any incident and reflected geometry in direct mode. Therefore, BRDF inversion was conveniently performed at custom source and view directions. Next, we obtained the black sky albedo (BSA) and white sky albedo (WSA) used for the albedo comparison in Section 3.3. These parameters are integrals over only the reflected hemisphere or bi-hemisphere and are given in Equations (14) and (15).

$$BSA({\theta }_i,\lambda )=\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\frac{\pi }{2}}R({\theta }_i,{\theta }_v,\phi ,\lambda )}}\sin {\theta }_v\cos {\theta }_{v}d{\theta }_{v}d\phi$$

(14)

$$WSA(\lambda )=2{\displaystyle {\int }_0^{\frac{\pi }{2}}BSA({\theta }_{{}_i},\lambda )}\sin {\theta }_{{}_i}\cos {\theta }_{{}_i}d{\theta }_{{}_i}$$

(15)

The angle index AFX, which is shown in Equation (16), was derived from kernel-driven models and has been applied to extract prior BRDF archetypes. Notably, $H_{\ker \_vol}$ and $H_{\ker \_geo}$ are the bi-hemisphere integrals of the two scattering kernels [38,39,40]. AFX was designed to effectively capture surface anisotropy, and values of AFX greater than 1.0 indicate a volume scattering-dominated bowl-like BRDF shape in the principal plane. In contrast, a dome-like BRDF shape associated with the main geometric-optical scattering will appear when the AFX is less than 1.0. The BRDF shape will be flat when the two scatterings are balanced. The index was introduced in this study to examine the potential scattering types of the PROSAIL model and their effects on model coupling.

$$AFX=\frac{WSA(\lambda )}{f_{iso}(\lambda )}=1+\frac{f_{vol}(\lambda )}{f_{iso}(\lambda )}\times H_{\ker \_vol}+\frac{f_{geo}(\lambda )}{f_{iso}(\lambda )}\times H_{\ker \_geo}$$

(16)

The RMSE and weight of determination ($WoD$) are generally used as evaluation criteria for kernel-driven models and are defined in Equations (17) and (18). In Equation (17), $R^{obs}$ and $R^{model}$ refer to $n$ observations and retrieved values. Roman et al. analyzed the retrieval accuracies of airborne and satellite observations at several spatial scales (30 m–500 m) for many land cover types and found that the average RMSEs of the RTLSR model were 0.01, 0.02, 0.05 and 0.04 in the 0.472-nm, 0.682-nm, 0.870-nm, and 1219-nm bands, respectively [41]. $WoD$ describes the amplification of noise from observations to inversion for angle sampling, as shown in Equation (18), where $K$ and $U$ are the kernel matrix and its weight vector, respectively [16,35]. The noise amplification of reflectance and albedo retrieval at given angles are included in $U$. Generally, three quality classes are established (0–0.75, 0.75–1.25 and >2), and a $WoD$ value less than one indicates less noise in the inversion process than in the original measurements.

$$RMSE=\sqrt{\frac{{\displaystyle {\sum }_{j=1}^n{(R^{obs}-R^{model})}^2}}{n-3}}$$

(17)

$$WoD={[U]}^T[{(K^{T}K)}^{-1}][U]$$

(18)

2.2. PROSAIL Canopy Reflectance Model

The PROSAIL model is a mature canopy radiation transfer model with the advantages of simplicity, robustness, and extensive validation. The latest 4SAIL canopy model considers hotspot effects, and the PROSPECT-5 model, which describes leaf optical reflectance and transmittance, is referenced in this study [6,7,14]. Additionally, simulation generation is a critical application for the PROSAIL model, and it can be in favor of the examination of the links among all the variables in the radiative transfer process [2].

At the leaf scale, the PROSPECT model can simulate directional-hemispherical reflectance and transmittance from 400 nm to 2500 nm for monocotyledonous, dicotyledonous and aging leaves [42]. Considering the rough surfaces of leaves, this model uses the maximum incident angle $\alpha$ to represent the diffused features. Meanwhile, an internal leaf is assumed to be composed of N parallel homogeneous layers. Reflectance and transmittance are calculated based on Equations (19) and (20), where ${\rho }_{90}$ and ${\tau }_{90}$ denote the reflectance and transmittance of every internal layer, and $R_{N-1,90}$ and $T_{N-1,90}$ are the total reflectance and transmittance of the internal N-1 layers, respectively. Finally, the PROSPECT model can be simplified to include only the plate transmission coefficient $\theta$ and the absorption coefficient $k$, which are given in Equations (21) and (22), where $k$ is related to the absorption $K_i$ of biochemistry component $i$, its weight $C_i$ and the absorption of albino leaf $k_e$ before 500 nm.

$${\rho }_l(\lambda )=R_{N,\alpha }={\rho }_{\alpha }+\frac{{\tau }_{\alpha }{\tau }_{90}{R}_{N-1,90}}{1-{\rho }_{90}{R}_{N-1,90}}$$

(19)

$${\tau }_l(\lambda )=T_{N,\alpha }=\frac{{\tau }_{\alpha }{T}_{N-1,90}}{1-{\rho }_{90}{R}_{N-1,90}}$$

(20)

$$\theta -(1-k)e^{-k}-k^2{\displaystyle {\int }_k^{\infty }{x}^{-1}{e}^{-x}dx=0}$$

(21)

$$k(\lambda )=k_e(\lambda )+{\displaystyle {\sum }_i\frac{C_i\cdot K_i(\lambda )}{N}}$$

(22)

The SAIL model is one of the earliest CR models based on the four-flux theory with 4 up- and down-scattering and absorption fluxes that are assumed for homogeneous, horizontal and infinitely extended vertical leaf layers [29]. The modified 4SAIL model is numerically robust and speed optimized, as given in Equation (23) [3].

$${\rho }_c(\lambda )=4SAIL(LAI,ALA,{\rho }_l(\lambda ),{\tau }_l(\lambda ),P_{soil}(\lambda ),SZA,VZA,\phi ,Hspot,SKYL),$$

(23)

where ${\rho }_c$ is the canopy reflectance, which couples ${\rho }_l$ and ${\tau }_l$ from the PROSPECT model. Other input parameters are shown in

Table 1. The input parameters of the PROSAIL model.

Parameter	Units	Common Value (range)
Leaf Scale
leaf structure parameter (Ns)	--	1.5 (1–3)
chlorophyll a and b contents (Cab)	μg/cm2	50 (0–100)
carotenoid content (Car)	μg/cm2	10 (5–30)
brown pigment content (Cbrown)	--	0
equivalent water thickness (Cw)	cm	0.015 (0.004–0.04)
leaf mass per unit leaf area (Cm)	g/cm2	0.009 (0.0019–0.0165)
Canopy Scale
LAI	m2/m2	3.5 (0–6)
average leaf angle (ALA)	degrees (°)	50 (0–90)
hot spot (Hspot)	mm	0.2 (0–1)
soil coefficient ($P_{soil}$)	--	0.1 (0–1)
diffuse/direct radiation (SKYL)	%	0 (0–100)
Observation geometry
SZA	degrees (°)	30 (0–90)
VZA	degrees (°)	0 (0–90)
relative azimuth angle ($\phi$)	degrees (°)	0 (0–180)

, and the reasonable ranges of parameters are based on previous studies of model intercomparison and sensitivity analysis [43,44]. Considering all the parameters, the extinction and scattering of the radiant flux of each leaf layer from the bottom to top of the canopy can be obtained, and CR can be calculated based on the cumulative energy transition.

2.3. Data and Experimental Design

Experimental design is the foundation of model comparison, and the overall design is shown in Figure 1. Almost all possible canopies simulated by the PROSAIL model can be considered via numerous simulations. In total, 20,000 vegetation parameter combinations were generated using the Saltelli periodic function [45] by the uniform sampling of nine leaf and canopy parameters over the ranges that are shown in Table 1 with reference to the experimental designs of previous studies [43,44], except for C_brown and SKYL. Considering that most leaves do not contain brown pigment [43] and the intrinsic BRDF is not necessarily related to the component of diffuse irradiance [46], the two remaining parameters were set to zero. Subsequently, the parameters were randomly combined to ensure diverse canopies. Then, we calculated the corresponding directional reflectances and albedos by direct simulation using the PROSAIL model in the red band (659 nm) and near-infrared (NIR) band (865 nm), where the angle distribution reflects the general geometry in remote sensing observations [47]. The intervals were 15°, 10° and 30° for the SZA (0–60°), VZA (0–80°) and $\phi$ (0–330°) with a total angle number of 397, respectively, based on typical angle ranges and experimental design [47,48]. POLDER can yield directional reflectance when SZA and VZA are close to 70° [49], and VZA ranges from nadir to 75° off-nadir for the airborne Cloud Absorption Radiometer (CAR) data set and MODIS satellite [33,41]. In addition, we set the maximum VZA to 80° for robust inversion ability, even at large angles close to 90° [17]. A $WoD$ of 0.0027 was observed, which indicated a good angle distribution with negligible noise, as discussed in Section 2.1.

As observations, the directional reflectances obtained from the PROSAIL model were then input into the kernel-driven RTCLSR/RTCLTR model to reconstruct the reflectance anisotropy. First, for the two adjustable hotspot factors C₁ and C₂ in the RTC model, the optimal values corresponding to the minimum RMSE of hotspot observations in each data set were determined [33]. Then, the optimal C₁ and C₂ factors were applied to the RTCLSR/RTCLTR model to retrieve the three kernel coefficients based on the database of 20,000 directional reflectances discussed in Section 2.1. There were 5390 and 1092 sets of data in the red and NIR bands, for which the volumetric scattering kernels were zero without hotspot effects. Subsequently, reflectances at the simulated angles and albedos (BSA at the five SZAs and WSA) were inverted based on the forward algorithm of the RTCLSR/RTCLTR model. Finally, we compared the directional reflectances and albedos of the two models based on the RMSE and bias. Notably, we called the comparison of directional reflectances “BRDF comparison” because the simulated data set of directional reflectance captured the main BRDF effects, which benefited from sufficient sampling at multiple angles. Before the comprehensive comparison, BRDF and albedo consistency was preliminarily examined for a typical canopy, and the BRDF effect was compared in the viewing hemisphere and principal plane. Additionally, we determined the reasons for model discrepancies by conducting statistical analyses based on vegetation parameters and RMSEs, including trend line and parameter distribution analyses based on whether error requirements were met, as noted in Section 2.1.

To further analyze the BRDF consistency of linking the two models, 21 sets of multiangle reflectances capturing a wide range of forest-grass-crop vegetation types were collected, where the LAI measurements were inferred from an LAI-2000 instrument. Most of these datasets are from PARABOLA radiometer measurements in the First International Satellite Land Surface Climatology Project (ISLSCP) Field Experiment (FIFE) [50], and the remaining data are from radiometer measurements by Kimes and Ranson [51,52,53]. The datasets were collected under nearly clear-sky conditions and have been widely used in the validation of the kernel-driven model [54,55,56] as well as extracting BRDF prior knowledge [40]. The central wavelength is 662 nm and 826 nm from PARABOLA radiometer from data 6 to data 18. For the remaining 8 sets of data that did not describe their central wavelength, we still used the 659 nm and 865 nm for PROSAIL simulation. There are 6801 BRF samplings altogether in which the VZA varies from 0 to 80°, and the SZA ranges from 16° to 82° with a total number of 126. Their detailed information is shown in

Table 2. The description of the 21 sets of BRDF field measurements.

Num	Land Cover	LAI (m2/m2)	Coverage (%)	Chlorosis (%)	Cloud (%)	SKYL (%)
1	Corn 1	0.65	25	0–5 4	0	0 6
2	Hardwood 2	4.20	79		0	0 6
3	Lawn grass 1	9.90	97	0–5 4	0	0 6
4	Orchard grass 3	1.10	50	30–40 4	0	0 6
5	Soybean 2	4.60	90	0–5 4	0	0 6
6	Mixed grass	1.50			10	15
7	Sand shinnery oak	1.75	60.2		5–10	(SZA, $\lambda$)
8	Grassland 3	1.39 ± 0.36		68.2 5	0	10
9	Grassland	1.55 ± 0.53		4.7 5	0	10
10	Grassland 3	1.19 ± 0.53		36.7 5	0	10
11	Dense grass-like vegetation 3	0.94 ± 0.30		65.9 5	0–5	10
12	Sparse grass 3	0.24 ± 0.05		42.9 5	0	10
13	Dense grass-like vegetation 3	0.60 ± 0.34		59.1 5	0–5	10
14	Grassland 3	1.91 ± 0.49		43.9 5	0–5	10
15	Sparse grass-like vegetation 3	0.12 ± 0.04		50.3 5	0	10
16	Grassland	3.59			5	8
17	Grassland	4.06			10	15
18	Cotton 2	4.00			0	10
19	Soybean 2	3.00 ± 0.50	72 ± 4		10–20	20
20	Soybean 2	3.90 ± 0.60	83 ± 3		1–20	15
21	Soybean 2	2.90 ± 0.40	99 ± 1		0	10

¹ The monocotyledons with N_s ranging from 1.0 to 1.5. ² The dicotyledons with N_s ranging from 1.5 to 2.5. ³ The senescent leaves with N_s larger than 2.5. ⁴ The percent of visual chlorosis. ⁵ The percent of dry weight of standing dead leaves among the total standing dry weight for both green and dead leaves. ⁶ The fraction of diffuse sky radiance had been corrected for these data.

, and the SKYLs for the remainder of the data were determined based on the observation of data 7 (SKYL varies with SZA and band), 16 and 17 (average SKYL), except for the first five groups which had been corrected based on the simulated diffuse fraction [53]. In addition, the N_s of some data were given based on vegetation classification and their fraction of dead leaves [57], and ${\rho }_l$ and ${\tau }_l$ for data 1 and 5 as well as ${\rho }_c$ for data 1 and 4 were observed. Because of the lack of vegetation measurements in these early campaigns, the other eight input parameters of the PROSAIL model were determined by looking for the optimal BRF fitting compared to the observations from the 20,000-parameter database mentioned above. Leaf inclination of the typical canopies for corn and soybean (data 1 and 5) presented a random distribution all over the range of $\phi$ [53], so we did not limit the scope of ALA during the parameter matching. Detailed information can be found in the footer of Table 2. Finally, we evaluated the two models based on their simulated or retrieved BRFs in comparison with the reflectance observations, and the resulting albedos were used for the cross-comparison.

3. Results

3.1. A Case Study for a Typical Canopy

In this section, we first provide an extreme example to parameterize a typical canopy as given in Table 1. The reflectance in the viewing hemisphere and principal plane that exhibits the most significant BRDF effect with SZA = 30° as well as the albedo, are shown in Figure 2. The polar plots have a common scale ranging from the minimum to the maximum for each row. The dome-shaped BRDF indicates that surface scattering is more obvious in the red band (Figure 2f), while volume scattering dominates in the NIR band, leading to a bowl-shaped BRDF (Figure 2i). The significant difference in scattering type is caused by high chlorophyll absorption in the red band and high leaf transmission in the NIR band, which lead to dominant single and multiple scattering effects, respectively. A slight overestimate can be seen in the forward direction, whereas backward scattering is slightly underestimated. This pattern is closely related to the solar principal plane because reflectance anisotropy is most significant in the principal plane, and discrepancies between models are easier to observe [31]. The hotspot-fits are somewhat better for the RTCLSR model than for the RTCLTR model, which demonstrates that the hotspot-corrected RossThick kernel (i.e., K_RTC) plays a better role in conjunction with the LSR kernel than the LTR kernel [30,36] in this case. The inconsistency can also be seen at large view angles in the red band, which shows that the RTCLTR model presents better fits than the RTCLSR model because the derivation of the LTR kernel includes the mutual shadowing effect [37,40]. These two models are highly consistent in the NIR band.

Figure 2d,e,g,h imply that bias between the kernel-driven model and PROSAIL models for the entire viewing hemisphere is less than 0.01 and 0.04 in the two bands. In addition, the absolute and relative RMSE of reflectances at 397 sampling angles between the two kinds of models are 0.0020 (1.10%, red, RTCLSR), 0.0017 (0.96%, red, RTCLTR), 0.0185 (2.53%, NIR, RTCLSR) and 0.0184 (2.53%, NIR, RTCLTR), respectively. Smaller RMSEs are found for the entire hemisphere with the statistic of 0.0012 (0.68%), 0.0010 (0.57%), 0.0082 (1.12%) and 0.0080 (1.10%). Notably, a significant difference can be seen at large relative azimuth angles especially for the NIR band, which may be related to the assumption of leaf distribution in the two models. There is an asymmetric BRDF shape along the principal plane due to the random leaf azimuth distribution for the PROSAIL model [29], whereas a symmetric shape is assumed in the kernel-driven model. In addition, a somewhat low-level accuracy at large angles for the semi-empirical Ross-Li model also contributes to such discrepancies [30,31]. Such a comparison for a typical example demonstrates that the BRFs by linking these two models are not necessarily consistent although the consistency between them is dominant.

The variability in BSA is similar to that of the BRFs, with a high value from the PROSAIL simulation when the SZA is greater than 40° and low values at smaller view angles (Figure 2). Notably, the largest deviations can arrive at 0.004 and 0.075 at 89° in the two bands. Despite these discrepancies, WSA exhibits good agreement between models, as shown in Figure 2. This result shows that once these two models are linked to evaluate the potentials in retrievals of biophysical/structural parameters using a simulation method, such a divergence between models should also be appropriately considered, although a generally good consistency exists between PROSAIL-generated BRFs and the fitting results by selecting an apt kernel-driven BRDF model.

3.2. General Consistency in BRFs between the Two Models

3.2.1. Comparison of Directional Reflectance Based on the Fit-RMSE and AFX Statistics

Using 20,000 sets of BRF simulations generated by the PROSAIL model that cover the main variations of the inputs from various biophysical/structural variables (Table 1), the comparison of BRFs between the RTCLSR and PROSAIL models is shown in Figure 3. The results of the RTCLTR model are generally similar to those of the RTCLSR model, so the results are not presented due to space limitations. Scatter plots are used to illustrate how the correlations are shown between the two models. The scatter plots are sliced into five equal-density levels in terms of the reflectance magnitude from 1 (maximum) to 5 (minimum) with the interval of 0.01 (Figure 3a,b). Figure 3a,b shows that the total R² is greater than 0.92, with negligible biases in the red and NIR bands, indicating that the two models are generally in good agreement for various canopies. 80% of the BRFs with density levels from 1 to 4 fall within the tolerance range although a few errors reach 0.11 and 0.17 in the two bands, as shown in Figure 3c and

Table 3. The AFX and fit-RMSE statistics between the RTCLSR/RTCLTR and PROSAIL models for diverse canopies.

Statistics	AFXRed	AFXNIR	RMSERed(LSR)	RMSENIR(LSR)	RMSERed(LTR)	RMSENIR(LTR)
Minimum	0.53	0.86	0.0000	0.0000	0.0000	0.0000
Maximum	1.56	1.60	0.1111	0.1745	0.1120	0.1745
Average	0.85	1.05	0.0092	0.0355	0.0103	0.0355

. Furthermore, 88.6% of the 20,000 fit-RMSEs are less than 0.02 in the red band, and 79.0% of canopies are within the tolerance range of 0.05 in the NIR band (Figure 3c). The total RMSEs are 0.0131 and 0.0436 (Figure 3a,b), with average fit-RMSEs of 0.0092 and 0.0355 in the two bands, respectively, as shown in Table 3.

The hemispherical distribution of average biases for the 20,000 sets of BRFs between the two models is calculated, where the angular bins refer to the BRDF sampling in Figure 1b and two of the SZAs that are shown in Figure 3d–i. When the SZA is within 30°, a slight underestimation of the PROSAIL-simulated hotspot BRFs by kernel-driven models can be found in the red band. As the SZA increases from 0° to 60°, model discrepancies tend to occur in large view zenith angles, and an obvious overestimation of the PROSAIL simulation BRFs by the kernel-driven BRDF model appears at the largest 80° VZA and 330° $\phi$ of the angle sampling. For the red band, the bias at large view zenith angles for the RTCLTR model is more significant than the result of the RTCLSR model (i.e., Figure 3d,e) except for the SZA of 60°, where fewer discrepancies are found for the RTCLTR model (Figure 3f,g). Unlike the red band, the presentation of model discrepancy between the PROSAIL model and the two kernel-driven models is generally similar in the NIR band. All results show that there is a general consistency in BRFs between these two models; however, few discrepancies in BRFs between two models mainly arise from hotspot direction and larger view zenith angles, somewhat depending on the variation of solar zenith angle. The RTCLTR model performs similarly in fitting these simulation BRFs generated by the PROSAIL model, although fitting results are somewhat better in a large view zenith angle (e.g., VZA > 65°). This result is generally consistent with the result examined in the case study in Section 3.1.

The AFX was applied to examine the scattering types among diverse canopies simulated from the PROSAIL model, as shown in Table 3 and Figure 4a. The AFX approximately satisfies the normal distribution clustered near 0.85 in the red band and the skewed distribution clustered at 0.95 in the NIR band (Figure 4a). To analyze the main BRDF properties in the two bands, 1.0 was taken as the critical value of AFX based on its characteristic [40]. Notably, 91.7% of AFXs are less than 1.0 in the red band, and the average value is 0.85, which implies that geometric-optical scattering is dominant. However, the numbers of data points are almost equal on either side of 1.0 and concentrated at approximately 0.95 when AFX is less than 1.0 in the NIR band. Additionally, the average value is 1.05. These results indicate primary isotropic scattering and volumetric scattering. Moreover, the range of AFX in the red band is wider than that in the NIR band, which is consistent with MODIS but over a narrower range [38,39,40]. The leading type of scattering in each band is related to high chlorophyll absorption in the red band and high leaf transmission in the NIR band, as discussed in Section 3.1. These factors strengthen single and multiple scattering. As a result, different scattering characteristics are observed for the two bands, and these characteristics indicate that the red band can capture more variability in the surface reflectance anisotropy than can the NIR band. The dominant volumetric scattering in the NIR band is theoretically caused by the multiple scattering of foliage leaves within canopies in the PROSAIL model. Furthermore, AFX variation has no significant impact on fit-RMSE in the red band (Figure 4b), while the RMSEs tend to increase as AFX increases in the NIR band (Figure 4c). This indicates that these models should be cautiously combined over canopies that represent extreme volume scattering effect (large AFX values), as is generally characterized by a dramatic BRDF shape, especially in the NIR band [40].

3.2.2. Model Discrepancies regarding Vegetation Parameters

To further explore the reason for the discrepancies of linking the RTCLSR and PROSAIL models, we examined the relationships between BRDF fit-RMSE of the two models and nine biophysical and background vegetation parameters. Figure 5a–f shows that the LAI and ALA can affect the RMSE in both the red and NIR bands and that $C_{ab}$ and $P_{soil}$ can affect the RMSE in the red band, all of which are sensitive parameters to CR [43,44]. Except for the sensitivity of fit-RMSE to the four input parameters of the PROSAIL model, other canopy structures exhibit subtle variations in RMSE throughout their ranges.

First, LAI is a critical factor that contributes to CR [23], and model discrepancies vary with LAI values. For a sparse or small-leaf canopy (LAI ≤ 1.0), RMSE exhibits an upward tendency, followed by a downward trend (LAI > 1.0) in the red band, as shown in Figure 5a, and a slight increase in the NIR band, as shown in Figure 5e. Meanwhile, large RMSEs are obtained when the LAI is less than 1.0 (Figure 5g). In contrast, as the LAI increases, the RMSE in the NIR band becomes more stable (Figure 5k). Notably, a relatively significant correlation can be found only for small LAIs in the NIR band, where R² is 0.2600. Consequently, model coupling is more applicable to low-density leaf or sparse canopies in the NIR band, whereas high-density leaf or dense canopies can be analyzed using the red band. These differences in relation to the LAI range may be associated with the different sensitivities of the LAI to CR. Notably, significant sensitivity can be observed at small LAIs in the two bands, while large LAIs are only sensitive in the NIR band [43,44].

Another essential canopy parameter is ALA, which is sensitive to CR in the NIR band [43,44] and displays similar tendencies in the two bands. Nevertheless, the fit-RMSE increases more significantly at large leaf angles in the NIR band, where R² is 0.5298, than in the red band (Figure 5b,f), and all the RMSEs that exceed 0.05 are associated with large ALAs (ALA > 60°) (Figure 5l). In addition, large ALAs also induce greater error in the red band (Figure 5b,h). Therefore, the kernel-driven model and PROSAIL model should be cautiously linked for vegetation estimation at large ALAs (mainly erectophile canopies), especially for the NIR band.

Contrary to the LAI, the fit-RMSE appears to decline at low-density leaf component $C_{ab}$ values and then levels off as $C_{ab}$ increases (Figure 5c). Simultaneously, the aggregation of large errors can be observed when $C_{ab}$ is less than 20 μg/cm² (Figure 5i). As a result, small $C_{ab}$ values are not a reasonable choice for the joint model.

The final parameter that contributes to the discrepancy between the kernel-driven model and PROSAIL model is $P_{soil}$. A subtle increase can be observed when the background soil becomes brighter, and consequently, most large errors accumulate when the soil exhibits high brightness. Therefore, in future studies, we should not focus on vegetation growing on bright soil, such as saline land, and should instead perform experiments on black soil.

3.3. General Consistency in Albedos between the Two Models

Land surface albedo is an essential factor in energy balance, as introduced in Section 2.1, and an absolute accuracy of 0.02~0.05 is usually required in many applications [58,59]. The GCOS report of 2016 (GCOS 2016) proposed the requirement that the absolute error of 0.0025 or relative deviation of 5% should be achieved within the next 10 years [60]. Additionally, land surface albedo is one of the main inversion parameters of the kernel-driven model, which is the integral of directional reflectance and can reflect the energy transfer. Therefore, we examined the consistency of linking the two models to provide some advice regarding albedo estimation. In this section, we do not show the results of the RTCLTR model because they exhibit little difference from those of the RTCLSR model. Despite the discrepancies in directional reflectance, Figure 6 shows good overall agreement between WSA and BSA at the 5 SZAs of simulations with R² values greater than 0.98 in the red and NIR bands. Moreover, less discrepancy is found in the WSA than in the BSA. Notably, a subtle overestimate of WSA is shown by the red band, whereas slight underestimates are found for the total BSA in both bands and for the WSA in the NIR band. Furthermore, the RMSEs of WSA between the two models are 0.0029 and 0.0156 in the two bands, which are within the relevant precision range. Although a wider range of AFX change is found in the red band, WSA and BSA exhibit better agreement in the red band than in the NIR band, which is consistent with the directional reflectance results in Figure 3a.

To assess the contribution of the SZA to the BSA in detail, the percentages of the albedo deviations in six ranges were calculated, as shown in

Table 4. The albedo deviation percentages between RTCLSR inversion and PROSAIL simulations in different accuracy ranges (%).

Band	Error	WSA	BSAtotal	BSA0 1	BSA15	BSA30	BSA45	BSA60
Red	[0, 0.02]	98.72	39.35	3.67	4.25	10.78	82.53	95.52
	[−0.02, 0)	1.26	59.30	93.25	94.57	89.21	17.33	2.13
	(0.02, 0.05]	0.01	0.60	0.51	0.29	0.00	0.00	2.21
	[−0.05, −0.02)	0.00	0.72	2.44	0.89	0.01	0.13	0.14
	>0.05	0.01	0.03	0.13	0.00	0.00	0.00	0.00
	<-0.05	0.00	0.00	0.00	0.00	0.00	0.01	0.00
NIR	[0, 0.02]	79.33	36.91	14.91	18.99	19.69	62.08	68.90
	[−0.02, 0)	18.66	49.60	64.32	63.99	79.60	24.23	15.82
	(0.02, 0.05]	0.81	4.77	9.33	13.14	0.22	0.17	1.00
	[−0.05, −0.02)	0.09	4.99	0.00	0.00	0.42	11.84	12.71
	>0.05	1.11	3.19	11.44	3.88	0.07	0.03	0.53
	<−0.05	0.00	0.54	0.00	0.00	0.00	1.65	1.04

¹ BSA₀ refers to the BSA when the SZA is zero (nadir), and BSA₁₅, BSA₃₀, BSA₄₅ and BSA₆₀ refer to the BSAs for SZAs of 15, 30, 45, and 60°, respectively.

and Figure 6c,f. Figure 6c,f clearly shows that most of the albedo deviations are less than 0.02 between these two models (red and blue areas), especially for the red band. However, some extreme errors ($\left|deviation\right|$ > 0.05) can be observed, particularly for the BSA when the SZA is 0° or 15° in the NIR band. The BSA varies greatly as the SZA changes, exhibiting the most similar error distribution to that of WSA at 60°, among the five SZAs tested. The transition from underestimation to overestimation can be clearly observed as the SZA increases, and this change is in agreement with the finding that the BSA exhibits an increasing trend as the SZA increases (Figure 2) [14]. This phenomenon is more obvious in the red band, and the largest difference is observed when the SZA is at nadir. Extremely high consistency is found at the turning point, at which almost all the canopies meet the 0.02 tolerance threshold.

Overall, the percentages of deviations less than 0.02 are 99.98% (WSA) and 98.65% (total BSA) in the red band and 97.99% and 86.51% in the NIR band. For a tolerance of 0.05, the percentages reach 99.99% and 99.97 in the red band, whereas 1.11% and 3.73%, respectively, of the canopies do not exhibit good agreement in the NIR band. In addition, the total relative deviations of WSA are 5.74% and 1.04% in the red and NIR bands, which can basically meet the requirement of GCOS 2016. We analyzed the reason given rise to these outliers especially for WSA in the NIR band (Figure 6d), where we found that the two sensitive parameters of LAI and ALA have a similar effect on WSA consistency compared to the result of reflectance in Figure 5e,f. In addition, little content of C_m seems to cause a larger difference. The good consistency in terms of the WSA may be explained by the compensatory effects of the BSA at different SZAs, which cause dramatic variance in the estimation accuracy as the SZA changes. Despite some extreme changes in the directional reflectance, the albedo is highly consistent, especially for the WSA, between the RTCLSR inversion and PROSAIL simulation results. Therefore, further research on albedo can be performed by coupling the two models and assessing their estimations.

3.4. Consistency in BRFs and Albedos between the Two Models Based on Field Observations

Similar to the comparison for the 20,000 PROSAIL simulation, the comparison of the simulated (PROSAIL) and retrieved (RTLSR) BRFs and albedos were also implemented based on 21 sets of BRF observations as shown in Figure 7, which is helpful in further exploring the BRDF consistency. For all the 6801 BRFs (Figure 7a–d), there is good agreement overall with the observation results. In addition, a better fitting of BRF observations for the RTLSR model can be seen than PROSAIL simulations in the red and NIR bands, which may be related to the uncertainty caused by vegetation parameter matching. Additionally, we preliminarily explored the relationship between angle index AFX and model error for each set of data (Figure 7e), from which we can see that the absolute biases between RTLSR retrievals and observations for all 21 data are close to the zero line with values less than 0.002, as well as in Figure 7f,g. Notably, bias presents a downward trend as AFX increases in the NIR band for PROSAIL simulations, which is inconsistent with the result of the 20,000 PROSAIL simulations (Figure 4c); therefore the influence of the AFX on model must still be investigated. In total, there are 14 sets of data where the absolute values of bias are less than 0.02 in the red band, and results for 12 sets of data are less than 0.05 in the NIR band for the PROSAIL simulations. For fit-RMSEs, the number is 19 and 14 in the two bands for RTLSR models, whereas smaller numbers are observed for the PROSAIL simulations, which are 7 and 5. The uncertainty from BRF measurements contribute to the model discrepancies, where each dataset was collected in a period by assuming a stable land surface [52]. The better fitting result from Ross-Li model is related to the extra error introduced by parameter matching of the PROSAIL model. In addition, the semi-empirical RTLSR model can be more flexible to fit the observations by the interpolation (or extrapolation) method than the physically-based PROSAIL model [56,61].

In terms of the relationship between the vegetation parameter and model discrepancy, similar patterns of fitting-error with the change of vegetation parameters are observed compared to Figure 5; this excludes soil reflectance, which seems to show no change. Bias between the PROSAIL simulation and observations is significant when LAI is less than 2.0 in the red band as shown in Figure 7f, which agrees with the results (Figure 5a). Nevertheless, the result for LAI is inconsistent in the NIR band (Figure 5e). In addition, bias is smaller when ALA is less than 30°, whereas larger bias is observed when ALA is greater than 50° in Figure 7f, which matches the early results well (Figure 5b,f). The statistics of average error for the four parameters mentioned in Figure 5 are shown in

Table 5. The statistics of the average error between the PROSAIL simulations/RTLSR inversions and observations for the four typical parameters in Figure 5, where the bias refers to the mean of their absolute values.

Error	LAI ≤ 1(>1)	ALA ≤ 40(>40)	Cab ≤ 20(>20)	Psoil ≤ 0.2(>0.2)
Bias_PROSAIL_Red	0.0185(0.0151) 1	0.0125(0.0167) 1	0.0183(0.0135) 1	0.0168(0.0149)
RMSE_PROSAIL_Red	0.0364(0.0258) 1	0.0190(0.0306) 1	0.0299(0.0277) 1	0.0249(0.0322) 1
RMSE_RTLSR_Red	0.0156(0.0096) 1	0.0110(0.0110) 1	0.0120(0.0085) 1	0.0103(0.0118) 1
Bias_PROSAIL_NIR	0.0463(0.0367)	0.0141(0.0448) 1
RMSE_PROSAIL_NIR	0.0743(0.0820) 1	0.0694(0.0827) 1
RMSE_RTLSR_NIR	0.0346(0.0455) 1	0.0570(0.0395)

¹ Error distribution is consistent with the results shown in Figure 5.

, where it is apparent that most of the error distribution is consistent with the results from the 20,000 sets of PROSAIL simulations in Figure 5. For a few of data that have a different trend with the results in Figure 5, these discrepancies may be due to the uncertainty mentioned above.

Finally, we inter-compared the albedo retrievals for the two models due to the lack of true values, including WSA and BSA at the 126 observed SZAs. As a whole, the values of WSA for RTLSR retrievals are lower than those of the PROSAIL simulations in the red band, whereas a better agreement is observed in the NIR band. For the fitting-error, all the R² values are greater than 0.76, and the total RMSE is 0.0212 and 0.0545 (with the relative deviation of 24.00% and 5.67%) for WSA in the two bands. In addition, there is a smaller RMSE for BSA. The albedo consistency for the 21 sets of observations is not as good as the result of the 20,000 sets of PROSAIL simulations, which may be caused by the uncertainty from BRF measurements and parameter matching as mentioned in the analysis for BRFs. In addition, the two models show a relatively good consistency without the effect of external uncertainty as shown in Figure 6, which confirms that these uncertainties are probably the major source that causes the somewhat large RMSE in Figure 8. In general, the two models have a generally good consistency in matching BRFs and albedos.

4. Discussion

In this study, the directional reflectance and albedo results of PROSAIL simulations as well as a variety of types of observations and their inversions using the hotspot-revised RTCLSR/RTCLTR model were compared for the first time, and good global agreement was observed, despite a few extreme discrepancies. The AFX was utilized to explore the scattering types among diverse canopies of PROSAIL simulations. Furthermore, the sensitivity of the model consistency to the vegetation structure was investigated to obtain a better understanding of the model linking. Additionally, the albedo, which is an essential factor in the models, was considered for further comparison to offer advice regarding albedo inversion by combining the two models.

The generally similar BRDF results of the RTCLSR model and the optimal RTCLTR model may be due to the robustness of the MODIS operational LSR model. Therefore, the commonly used RTCLSR model can be used in future model coupling. The case study results involving a common canopy exhibited good overall agreement, indicating that model coupling can be performed when only typical vegetation is considered. Despite the extreme error in the directional reflectance (Figure 3), 88.6% and 79.0% of canopies can meet tolerances of 0.02 and 0.05 in the red and NIR bands, respectively. On one hand, the kernel-driven model has proven to be reliable for various canopies [56,62,63], and therefore such a consistency is not surprising. On the other hand, however, discrepancies between these two models remain and may be caused by the differences in model theory, as discussed in the introduction. In addition, the diverse combinations of vegetation parameters required to cover complete canopies may be responsible for these differences because unreasonable combinations that do not exist in the real world may be included. However, no better way of performing a complete analysis has been developed to date, and the method of random combination has been widely applied in comprehensive analyses [43,44]. Therefore, further studies of the interactions among vegetation variables are necessary to guide the experimental design. Nonetheless, we have obtained an overall understanding of the consistency between the two models within the capacity of the PROSAIL model, and these findings will facilitate further studies on joining these two models. In terms of the field measurements, the RTLSR model presents better fitting compared to the PROSAIL model, which may imply that the semi-empirical kernel-driven models have more potentials in flexibly fitting actual multiangle observations than the physical models that require the parameterization of various biophysical/structural and spectral variables as inputs. In addition, more discrepancies were observed in the result for the field data than that of the 20,000 simulations, which may be related to the uncertainties during the BRF measurements and vegetation parameter matching.

To further explore the reasons underlying these discrepancies, the links between the vegetation structure and model differences were analyzed in detail, and canopy parameters, including LAI, ALA and $P_{soil}$, were shown to be sensitive in the red band. Similarly, LAI and ALA can also affect the fit-RMSE in the NIR band. At the leaf level, variations can only be observed as $C_{ab}$ changes. This result suggests that the largest effect on model discrepancy occurs at the canopy level, especially for ALA, and that the canopy structure can significantly affect reflectance in the red and NIR bands, as shown in previous studies [1]. Because of the equal variation derived from the uniform sampling of every parameter, the sensitivity analysis can be considered reliable. Hence, the analyses of canopies with small LAI and C_ab values and those with large ALA and $P_{soil}$ values suggest that we should be cautious when combining the two models in the red band and that further examination is required for both large LAIs and erectophile canopies in the NIR band. Notably, a real vertical leaf is limited [29]. This limitation may affect the validation of PROSAIL model simulations and, subsequently, comparisons such as those performed in this study. Moreover, large ALAs (erectophile canopy) generally cause a dramatic change in the BRDF shape [29], for which a large hotspot is generally observed. Our group is working on a hotspot correction for the geometric-optical scattering kernel, and a reexamination of these discrepancies is likely necessary using the updated model.

Parameter inversion is a potential application of linking these two models. The two models used in this study benefit from both the complete description of the vegetation structure in the PROSAIL model and the robust BRDF simulation and inversion of the kernel-driven model. In this study, we examined the parameter consistency in terms of albedo as an example. The highly consistent albedo reflects reliable model coupling for estimation and evaluation. Limited to the five SZAs used for BSA comparison, the BSA consistency must be further examined for other SZAs. Similarly, good consistency was observed for the 21 sets of field BRDF observations, in which the lower agreement may be due to error that occurs during the post-processing of the vegetation parameter matching and the uncertainty of from BRF measurements themselves. Based on the parameter matching, the retrieved and measured LAIs were compared. The result shows good consistency with a high R² of 0.8584 and a total RMSE of 0.95, which can probably meet the precision of current LAI products (around ±1.0) from MODIS and CYCLOPES sensors [64]. Such an attempt for the LAI retrieval may initially demonstrates the feasibility of PROSAIL inversion by linking the two models. The three coefficients $f_{iso}$, $f_{vol}$ and $f_{geo}$ included in the MODIS BRDF product are promising for estimating vegetation parameters, such as in the CI [25,26] and canopy height [28], as discussed in the introduction. In addition, the AFX has been applied to extract prior BRDF archetypes and can capture the main BRDF variety, which is potentially related to the vegetation architecture [38,39,40]. Thus, these BRDF factors could be added as constraints for parameter retrieval using the PROSAIL model in the future, where the relationships between these BRDF parameters and the vegetation structure can be further studied. Notably, the model consistency analysis in this study provides meaningful guidance for such studies. We preliminarily analyzed the sensitivity of retrieved Ross-Li model coefficients against PROSAIL parameters by using the Extended Fourier amplitude sensitivity test (EFAST) method. The result is consistent with Figure 5, which suggests that sensitive parameters on BRDF effect occurs at the canopy level, especially for the ALA. Additionally, real data, such as sufficient measurements of LAI and ALA as well as other vegetation parameters, can be used to verify the results of these comparisons during their estimations. The strategy of comparison of BRDF for linking different models proposed in this study is still applicable for many other physical BRDF models beyond the PROSAIL model, and it would be promising to estimate various vegetation structure related to the physical models from satellite BRDF data by considering the BRDF consistency of model linking.

5. Conclusions

In this study, the potential in linking two widely used models (i.e., a semi-empirical, hotspot-revised, kernel-driven model, including the RTCLSR and RTCLTR models, and the physical PROSAIL model) was comprehensively investigated in terms of BRDF and albedo. A typical canopy was first exemplified in the comparison in viewing hemispherical and principal planes to examine the potential in linking the two models. In total, 20,000 sets of parameter combinations were applied in a comprehensive analysis based on the PROSAIL model and were used to generate the PROSAIL multiangle reflectance and albedo. Then, BRDF and albedo values were obtained using the kernel-driven model. Finally, the reflectance and albedo values of the two models were compared, and the correlations between vegetation parameters and model discrepancy were further explored. To further strengthen the investigation, 21 sets of field BRDF observations that cover a variety of land surface classes were also introduced into the study.

Overall, similar results were observed for the RTCLTR and RTCLSR models. This finding indicates that the MODIS operational LSR model provides a robust inversion technique, although the LTR model is considered to be the optimal version of the LSR model [36]. Globally, the two models agree well in terms of directional reflectance for diverse canopies (R² > 0.92), where 88.6% and 79.0% of fit-RMSEs are within tolerances of 0.02 and 0.05, respectively, in the two bands. This indicates great potential for linking these two models for potential applications, such as for developing inversion procedures to accurately retrieve vegetation properties from the satellite BRDF parameter product generated by the kernel-driven BRDF model. Limitations may arise from the few inconsistencies in BRFs between the two models in a large VZA beyond 70°, which somewhat depends on the variation of SZAs. In terms of reflectance anisotropy, geometric-optical scattering was dominant in the red band, whereas volumetric scattering was usually dominant in the NIR band. Model differences slightly increase for extreme case of volume scattering (large AFX values), especially in the NIR band. Additionally, attention should be paid to small LAI and C_ab values and large ALA and $P_{soil}$ values in the red band, and large LAI and ALA values in the NIR band, which tend to result in more errors when linking the two models, especially for erectophile canopies. Great potential in linking the two models should arise based on the good agreement in albedo (R² > 0.98), and 99.98% and 97.99% of the WSA deviations were less than 0.02 in the two bands. At an SZA of 30°, the BSA tends to perform best among all simulations. In addition, a similar result was found for the 21 sets of field BRDF observations, which presents more consistency between the two models than the result of the 20,000 PROSAIL simulations. This may contribute to the real measurements of BRDF and vegetation structure.

The overall good consistency between the two models suggests a great potential for linking these two models for parameter estimation, especially from the BRDF parameter product generated by the kernel-driven BRDF model, except for extreme cases that involve a few types of canopies. Therefore, future research should focus on the application of prior BRDF knowledge derived from the kernel-driven model, such as global BRDF parameters (e.g., MCD43A) and AFX BRDF archetypes, in conjunction with the PROSAIL inversion to devise an inversion relationship between vegetation structure and model BRDF parameters, based on the findings in this study.