Limitations and Improvements of the Leaf Optical Properties Model Leaf Incorporating Biochemistry Exhibiting Reflectance and Transmittance Yields (LIBERTY)

Abstract

Leaf Incorporating Biochemistry Exhibiting Reflectance and Transmittance Yields (LIBERTY) models the effects of leaf biochemical concentrations on reflectance spectra on the basis of Melamed theory, which has several limitations. These are: (1) the radiation components are not treated satisfactorily; (2) the directional changes of both particle and sublayer scattering ratios are not considered; and (3) the boundary constraint which makes needle leaves different from broadleaves is not included. Proofs of these limitations as well as theoretical improvements are given in this study. Global sensitivity analysis (SA) of three models: the original LIBERTY, our improved LIBERTY (LIBERTY_im) and The optical PROperties SPECTra model (PROSPECT) suggests that compared with LIBERTY, the global reflectance and transmittance of LIBERTY_im are more sensitive to diametrical absorbance $\alpha d$—a parameter related to leaf biochemistry. Moreover, the global reflectance and transmittance of LIBERTY_im and PROSPECT had similar sensitivity patterns to the input variables, demonstrating indirectly the validity of our improvements over LIBERTY. However, neither LIBERTY nor LIBERTY_im considers boundary constraints, which limits their applications in modelling needle leaf optical properties. We introduced a particle string model, which might be used to simulate needle leaf optical properties in the future.

1. Introduction

As the main light-harvesting organ of vegetation, plant leaves undoubtedly play a vital role in the carbon dynamics of terrestrial ecosystems. Great interest has been shown by botanists, spectroscopists and remote sensing scientists and practitioners to accurately model leaf optical properties. The main difficulties, as pointed out by Wang, et al. [1], come from the intricate underlying structures and their complicated interactions with lights. The approaches for leaf diffuse optical properties can be classified into six categories: (1) plate models; (2) compact spherical particle models; (3) N-flux models; (4) radiative transfer equation; (5) stochastic approach; (6) ray tracing models. These models will not be discussed here since they have been fully reviewed by Jacquemoud and Ustin [2]. We will only focus on the spherical particle model Leaf Incorporating Biochemistry Exhibiting Reflectance and Transmittance Yields (LIBERTY), which was developed to model the effects of leaf biochemical concentrations on reflectance spectra for needle leaves particularly.

LIBERTY was proposed by Dawson, Curran and Plummer [3], and originates from Melamed theory [4] and Benford theory [5]. As illustrated in Figure 1, it treats a leaf cell as a homogenous particle and a leaf is regarded as an ensemble of particle layers. More specifically, Melamed theory is used to calculate the reflectance of infinite particle medium $R_{\infty }$ and the sublayer reflectance $R_s$. Then, the sublayer transmittance $T_s$ can be obtained with $R_{\infty }$ and $R_s$. Finally, the number of sublayers $N$ together with $R_s$ and $T_s$ are input into Benford theory to calculate the global reflectance $R_N$ and transmittance $T_N$. However, Moorthy, et al. [6] found the LIBERTY model was ineffective in retrieving needle-level chlorophyll concentration. The reason is that several limitations exist when applying Melamed theory to model leaf optical properties of needle leaves. The limitations are as follows: (1) the radiation components are not treated satisfactorily; (2) the directional changes of both particle and sublayer scattering ratios are not considered; and (3) the boundary constraint which makes needle leaves different from broadleaves is not included.

This paper is structured as follows. We first demonstrate the existence of limitations in LIBERTY in Section 2. Section 3 discusses these limitations in detail. Theoretically-corrected Melamed theory is given in Section 4. In this section, we also compared LIBERTY and The optical PROperties SPECTra model (PROSPECT), and the reflectance and transmittance of a particle string is worked out, which we suggest may be able to model leaf optical properties of needle leaves in the future. Section 5 describes data ranges and sensitivity analysis method. Results are shown in next section. Section 7 discusses in detail the differences between three models and the potentials of the particle string model. Conclusions are drawn in Section 8. The study is not a specific critique. The aim of this study is to point out the common obstacles met in the way of needle leaf optical properties modelling, which may be helpful for future studies.

2. Existence of Limitations

2.1. The Reflectance of Infinite Particle Medium ($R_{\infty }$) Is Inconsistent with the Standard Form

The equation of $R_{\infty }$ in the original LIBERTY model is:

$$R_{\infty }=\frac{2x{m}_e+x\left(1-2x{m}_e\right)\tau \left(1-m_e{R}_{\infty }\right)}{\left(1-m_e{R}_{\infty }\right)-\left(1-x\right)\left(1-m_e\right)\tau R_{\infty }}$$

(1)

where $x$ is the particle backscattering ratio, $m_e$ is the external reflectance, $\tau$ is the transmittance of a particle, $R_{\infty }$ is the reflectance of infinite particle medium. Their equations can be found in Nomenclature.

However, typographical and other errors in the original physical formulation of the theory have been found and a corrected version given by Mandelis, et al. [7] is:

$$R_{\infty }=2x{m}_e+x\left(1-2x{m}_e\right)\tau \left[\frac{1-m_e{R}_{\infty }}{\left(1-m_e{R}_{\infty }\right)-\left(1-x\right)\left(1-m_e\right)\tau R_{\infty }}\right]$$

(2)

In LIBERTY, the infinite reflectance $R_{\infty }$ and the sublayer reflectance $R_s$ are used to calculate the sublayer transmittance $T_s$. Then, the number of sublayers $N$, the sublayer reflectance $R_s$ and transmittance $T_s$ were input into a multiple plate model to obtain the global reflectance and transmittance. The infinite reflectance $R_{\infty }$ is the limiting case of a multiple plate model when the number of sublayers is infinite. Therefore, $R_{\infty }$ must be consistent with the standard form otherwise the accuracy of the calculated $T_s$ will be compromised. The reflectance $R_{N+1}$ and transmittance $T_{N+1}$ of a multiple plate model with $N+1$ sublayers can be expressed as the following equations:

$$R_{N+1}=\frac{T_s^2{R}_N}{1-R_s{R}_N}+R_s$$

(3)

$$T_{N+1}=\frac{T_s{T}_N}{1-R_s{R}_N}$$

(4)

Both Benford theory [5] used by LIBERTY, and the Stokes theory [8] used by PROSPECT are built on the two equations above, thus they are essentially the same. From Equation (3), we can get the standard form of $R_{\infty }$:

$$R_{\infty }=\frac{T_s^2{R}_{\infty }}{1-R_s{R}_{\infty }}+R_s$$

(5)

Dawson, Curran and Plummer [3] assume the sublayer reflectance ($R_s$ i.e., rreflectance of elementary particle layer) is:

$$R_s=2x{m}_e+x\left(1-2x{m}_e\right)\tau$$

(6)

where $x$ denotes the sublayer backscattering ratio, $m_e$ is the external reflectance coefficient and $\tau$ is the particle transmittance. However, Equation (2) cannot be converted to the form of Equation (5) after taking Equation (6) into it.

In PROSPECT, the treatments of incident and scattered radiation are different because the maximum incident angle of incident radiation is assumed to be smaller than 90°, but for the scattered radiation the maximum incident angle is 90°. LIBERTY may be built on the same philosophy since it considers the sublayer scatterings on the particle surface for incident radiation only when calculating the sublayer reflectance and transmittance. In this case, Equation (3) would be changed into:

$$R_{N+1}=\frac{T_1{T}_s{R}_N}{1-R_s{R}_N}+R_1$$

(7)

where $R_1$ and $T_1$ denote the sublayer reflectance and transmittance for incident radiation, and $R_s$ and $T_s$ are sublayer reflectance and transmittance for the scattered radiation, respectively. The form of $R_{\infty }$ is:

$$R_{\infty }=\frac{T_1{T}_s{R}_{\infty }}{1-R_s{R}_{\infty }}+R_1$$

(8)

If the treatment of incident radiation is assumed to be unique, then $R_1=2x{m}_e+x\left(1-2x{m}_e\right)\tau$, $T_1=\left(1-x\right)\left(1-2x{m}_e\right)\tau$, $R_s=m_e+x\left(1-m_e\right)\tau$, $T_s=\left(1-x\right)\left(1-m_e\right)\tau$. After taking $R_1$, $T_1$, $R_s$ and $T_s$ into Equation (8), we cannot get the form like Equation (2) either.

Mandelis, Boroumand and Bergh [7] physically corrected Melamed theory and gave a different equation for $R_{\infty }$

$$R_{\infty }=2x{m}_e+x\left(1-2x{m}_e\right)\tau \left[\frac{\left(1-x\right)\left(1-m_e\right)\tau R_{\infty }+1-m_e{R}_{\infty }}{\left(1-m_e{R}_{\infty }\right)-\left(1-x\right)\left(1-m_e\right)\tau R_{\infty }}\right]$$

(9)

Unfortunately, it cannot be converted to the standard form either.

2.2. The Sublayer Transmittance $T_s$ Is Not Equal to the Theoretical Value

As mentioned in Section 1 and Section 2.1, the infinite reflectance $R_{\infty }$ together with the sublayer reflectance $R_s$ is used to calculate the sublayer transmittance $T_s$. In LIBERTY, $T_s$ is calculated using Equation (5), that is:

$$T_s={\left[\frac{\left(R_s-R_{\infty }\right)\left(1-R_{\infty }{R}_s\right)}{R_{\infty }}\right]}^{1/2}$$

(10)

If there is no underlying leaf material whose backscattering will contribute to the total reflective radiation ($R$ in Figure 1), theoretically, the sublayer transmittance will be Component 3 in Figure 1, that is:

$$T_s=\left(1-x\right)\left(1-2x{m}_e\right)\tau$$

(11)

The sublayer transmittances calculated using two equations above are supposed to be equal to each other. Since $x$, $m_e$ and $\tau$ are determined by four parameters (see Nomenclature): the relative refractive index $n$, the particle backscattering ratio $x_u$, the particle forward scattering ratio $x_d$ and the diametrical absorbance $\alpha d$, both $T_s$ in Equations (10) and (11) can be represented as a function of these four parameters i.e., $T_s\left(n,\alpha d,x_u,x_d\right)$. We set $n, x_u,x_d$ to fixed values of 1.45, 0.2, 0.6 respectively, thus $T_s\left(n,\alpha d,x_u,x_d\right)$ becomes a function of the diametrical absorbance $\alpha d$. As shown in Figure 2, the sublayer transmittance of original LIBERTY is underestimated by contrast with theoretical values. This discrepancy can be attributed to the invalid equation of $R_{\infty }$.

2.3. The Fraction of Incident Radiation Entering the Sublayer is Overestimated

In LIBERTY, the fraction of incident radiation entered into the first sublayer is assumed to be $1-2x{m}_e$ (Component 1 in Figure 1); we can change the form into:

$$1-2x{m}_e=1-m_e+m_e-2x{m}_e=1-m_e+\left(1-2x\right)m_e$$

(12)

where $1-2x$ denotes the sublayer attenuation ratio since the sublayer backscattering ratio is equal to the forward scattering ratio. The attenuated part has been added to the radiation penetrated into the particle. Even if Component 1 is correct, the following reflections on the lateral particles (Components 6, 9, 12) will be the same. However, when the reflected radiation by underlying particle bulk strikes on the lateral particles, the reflectance turns into $m_e$. As mentioned in Section 2.1, LIBERTY may have special treatment for incident radiation since it is collimated rather than scattered radiation. However, this cannot explain the overestimated radiation component in the equation above.

3. Limitations of LIBERTY

3.1. The Directional Changes of Both Particle and Sublayer Scattering Ratios Were Not Considered

One of the limitations that exists in LIBERTY is that the directional changes of the scattering ratios in both particle and sublayer scales are not considered. If the particle is homogeneous, its backward, lateral and forward scattering ratios ($x_u$, $x_a$ and $x_d$) will be vectors whose directions depend on that of incident radiation. Figure 3 gives a more intuitive representation. Since the sublayer backscattering ratio $x$ is related to $x_u$, $x_a$ and $x_d$ (see Nomenclature), the direction of $x$ will be determined by that of incident radiation as well. However, in LIBERTY the direction of $x$ is fixed upward. As noted in Figure 1 in reddish tone, $x$ is used in Component 2, 5, 8, 11. However, Component 1 is different from Components 5, 8, and 11 because the radiation is incident from above while for Components 5, 8, and 11 the radiation is from below.

3.2. The Radiation Components Were Not Treated Satisfactorily

Several radiation components in Figure 1 need proper treatment. The first is Component 1. As mentioned in Section 2.3, the original Melamed theory assumes that there is $2x{m}_e$ of incident radiation reflected on particle surface; the remaining part $1-2x{m}_e$ will propagate into the particle. First let us review the definition of $x$. It represents the backward fraction of outgoing radiation from a sublayer. Figure 4 shows a more intuitive representation. Radiation is incident on the left particle (the entrance particle) from above, and the radiation outgoing from this particle is assumed to be a unit, out of which $x_a$ will propagate to lateral particles, and repeat the same transfer process. Reflections on the adjacent particle surface are not considered, because in real scenarios cells are connected with each other by cell walls, so intercellular radiation will penetrate entirely into cells by multiple reflections between the cell walls. Apparently, $x$ takes account of influences of the lateral particles in the same sublayer. The final equation of $x$ is:

$$x= x_u+x_a\tau x_u+{\left(x_a\tau \right)}^2{x}_u+\dots =x_u/\left(1-x_a\tau \right)$$

(13)

Melamed [4] assumes $x_u\approx x_d$ at low absorption coefficients where $\alpha d<1$. Hence, the sublayer backscattering and forward scattering ratios are equal to each other. There will be $2x{m}_e$ of incident radiation reflected on the particle surface because “the reflected component $m_e$ of the external incident ray is scattered over an angle $2\mathsf{\pi }$ steradians therefore $x_d$ is zero and $x_u$ is twice the radiation which emerges from within a cell” [3]. However, the reason does not hold. As shown in Figure 4, for the middle particle which receives $x_a$ of outgoing radiation from the left particle, the scattering still happens over an angle of $4\mathsf{\pi }$ steradians. In addition, since $x$ has considered the shadowing of the lateral particles, as explained in Equation (12), the attenuated part has been added into the radiation penetrating into the particle. Obviously, $2x{m}_e$ does not conform to the reality. Then, how much incident radiation enters the particle on which it strikes? We believe there is no need to consider the scattering and shadowing of the lateral particles when calculating the penetrated part of incident radiation. Therefore, there will be just $m_e$ of incident radiation reflected on the surface and $1-m_e$ entering into the particle.

The second radiation component is Component 3. Dawson, Curran and Plummer [3] assume that the sublayer backscattering ratio is $x$, therefore there will be $1-x$ of radiation scattered forward. However, according to the definition of $x$, the shadowing of the lateral particles has been considered, which means a fraction of the total radiation will be attenuated inside the sublayer during radiation transfer. In Figure 4, the radiation transferred into the middle particle is $x_a$, out of which $1-\tau$ will be attenuated in this particle ($\tau$ is the transmittance of a particle). The attenuation will also happen in the right particle, and $1-x$ has included this attenuated radiation. If the particle forward scattering ratio $x_d$ is equal to $x_u$, then there will be $x$ of radiation scattered forward by the sublayer.

3.3. The Morphology of Needle Leaves Was Not Included

LIBERTY is developed to model optical properties of needle leaves but it does not include the morphological structures of needle leaves which makes them different from broadleaves. Needles are narrow and thick, and needle structures with adaxial and abaxial surfaces are neither parallel nor necessarily flat planes [9]. Figure 5 shows an anatomical structure of a typical needle. The limited boundary of single particle sublayer has not been considered in LIBERTY. The number of particles in a sublayer is assumed to be infinite in LIBERTY. Actually, if the boundary constraint is not included, the main difference between LIBERTY and PROSPECT is the sublayer morphology. The sublayer in LIBERTY is composed of compact particles while in PROSPECT it is a flat plate. The two global reflectance calculation methods—Benford theory and Stokes theory, used by LIBERTY and PROSPECT respectively—are essentially the same.

4. Improved Version of LIBERTY

Based on the issues raised above associated with LIBERTY, this section will propose an improved version of LIBERTY, which is called LIBERTY_im for convenience.

4.1. Reassessment of the Sublayer Backscattering and Forward Scattering Ratios

In LIBERTY, the particle backscattering ratio $x_u$ is assumed to be equal to the particle forward scattering ratio $x_d$. However, the prerequisite of this assumption, that is $\alpha d<1$, does not always hold (please refer to

Table 1. Ranges of input parameters. $n$ is the refractive index, $\alpha d$ is the diametrical absorbance, $N$ is the number of sublayers, $x_u$ is the particle backscattering ratio.

Input Parameters	Minimum Value	Maximum Value
$n$	1.2708	1.5295
$\alpha d$	0.0021	8.4985
$N$	1.0927	3
$x_u$	0	0.5

). Hence, the $x_u$ and $x_d$ are not necessarily equal to each other in LIBERTY. As a result, the particle forward scattering ratio $x_d$ must be input into the model as an additional parameter and the sublayer backscattering and forward scattering ratios are not necessarily equal to each other either. The sublayer backscattering and forward scattering ratios are denoted as $X_u$ and $X_d$ in LIBERTY_im.

As explained in Section 3.1, it is more reasonable to consider the directional changes of the sublayer backscattering and forward scattering ratios (as shown in Figure 3). When the directions of the particle scattering ratios (i.e., $x_u$, $x_a$ and $x_d$) are taken into account, the calculations of the sublayer scattering ratios are similar to Figure 6. Radiation enters the sublayer through the left particle, which is called the entrance particle. For the middle particle, the incident radiation comes from the left, that is $x_a$. When the directions of the particle scattering ratios change, the proportions of scattered radiation received by the lateral particles change as well since the receptors for $x_u$ and $x_d$ are sublayers (see Figure A1a,b). We define two direction change coefficients $p$ and $p_a$ to represent the proportions of $x_u$($x_d$) and $x_a$ received by lateral particles when the directions of particle scattering ratios change. The laterally scattered radiations (the middle and right particles in Figure 6) are assumed to undergo a similar transfer process, so they have been put together. The reassessed value for the sublayer backscattering ratio $X_u$ can be derived by an infinite series summation that converges to:

$$X_u=x_u+\frac{\left[1-p+\left(p-p_a\right)x_a\right]x_a\tau }{2-2\tau \left[p-\left(p-p_a\right)x_a\right]}$$

(14)

In the equation above, a transformation is used: $x_u+x_d=1-x_a$. Similarly, the sublayer forward scattering ratio $X_d$ can be calculated by following equation:

$$X_d=x_d+\frac{\left[1-p+\left(p-p_a\right)x_a\right]x_a\tau }{2-2\tau \left[p-\left(p-p_a\right)x_a\right]}$$

(15)

Then $x_u$ to $X_u$ and $x_d$ to $X_d$ ratios denote the contributions of the entrance particle to the sublayer backscattering and forward scattering ratios, respectively. When considering the direction changes of particle scattering ratios, $x_u$ and $x_d$ are exactly the first left items on the right-hand side of Equations (14) and (15). Thus, the second left items are same and represent the contributions of the lateral particles to sublayer backscattering and forward scattering ratios.

4.2. Reassessment of the Sublayer Reflectance and Transmittance

The sublayer backscattering ratio $X_u$ and forward scattering ratio $X_d$ do not include the reflection and transmission on the surface of the entrance particle (i.e., rthe left particle in Figure 6). After we get the formulae of $X_u$ and $X_d$, the equation of sublayer reflectance and transmittance can be derived by:

$$R_s=m_e+\left(1-m_e\right)\tau X_u$$

(16)

$$T_s=\left(1-m_e\right)\tau X_d$$

(17)

Stokes theory instead of Benford theory is used to calculate the global reflectance and transmittance. The philosophy of Benford theory is to break down thickness into fractional part ranging from 1 to 2 and the remaining integer part. Actually, the difference in reflectance calculated from these two theories is so tiny that it can be neglected (results not shown). This is due to the same premise on which these two theories are built as mentioned in Section 2.1. Compared with Benford theory, Stokes theory is simpler. The global reflectance $R_N$ and transmittance $T_N$ of $N$ homogeneous layers are:

$$\frac{T_N}{a-a^{-1}}=\frac{R_N}{b^N-b^{-N}}=\frac{1}{a{b}^N-a^{-1}{b}^{-N}}$$

(18)

where, $a=\left(1+R_s^2-T_s^2+\Delta \right)/2T_s$, $b=\left(1-R_s^2+T_s^2+\Delta \right)/2T_s$, $\Delta =\sqrt{\left(1+R_s+T_s\right)\left(1-R_s+T_s\right)\left(1+R_s-T_s\right)\left(1-R_s-T_s\right)}$.

4.3. The Reflectance and Transmittance of a Particle String

It is difficult to consider boundary constraints for a sublayer. We just give the solution of a particle string.

Suppose a limited number (M+N) of particles are arranged compactly in a straight line and unit light is incident isotropically from the left of the system (shown in Figure 7a). Since the system is composed of particles with the same and homogeneous quality, we can regard the M particles as a block and remaining N particles as another block (shown in Figure 7b). If the particle scattering ratios (i.e., $x_u$, $x_d$ and $x_a$) are constants, $x_u$ amd $x_d$ can be treated as the sublayer reflectance and transmittance, thus the global backscattering ratio (U) and forward scattering ratio (D) of the particle string can be calculated by Equation (18). Specific equations are given in Nomenclature.

In order to get the equation of lateral scattering ratio, we list three directional fractions of the nth transfer in

Table 2. Statistics of three directional fractions of outgoing radiation in the nth radiation transfer in Figure 7b. $U$ and $D$ are the backscattering and forward scattering ratios of the particle string, respectively. $R$ and $T$ are the reflectance and transmittance of the particle string, respectively.

Count	Backward *	Forward *	Reflected *	Next **
1	$U_M$	$D_M$	$R_M$	$D_M$
2	$D_M{U}_N$	$D_M{D}_N$	$D_M{R}_N$	$D_M{U}_N$
3	$D_M^2{U}_N$	$D_M{U}_N{U}_M$	$D_M{U}_N{T}_M$	$D_M{U}_N{U}_M$
4	$D_M{U}_N^2{U}_M$	$D_M{U}_N{U}_M{D}_N$	$D_M{U}_N{U}_M{R}_N$	$D_M{U}_N^2{U}_M$
5	$D_M^2{U}_N^2{U}_M$	$D_M{U}_N^2{U}_M^2$	$D_M{U}_N^2{U}_M{T}_M$	$D_M{U}_N^2{U}_M^2$
6	$D_M{U}_N^3{U}_M^2$	$D_M{U}_N^2{U}_M^2{D}_N$	$D_M{U}_N^2{U}_M^2{R}_N$	$D_M{U}_N^3{U}_M^2$
7	$D_M^2{U}_N^3{U}_M^2$	$D_M{U}_N^3{U}_M^3$	$D_M{U}_N^3{U}_M^2{T}_M$	$D_M{U}_N^3{U}_M^3$

* directions are relative to incident radiation. ** radiations left for next transfer.

. The total laterally-scattered radiation $R_{M+N}$ equals the sum of the “Next” column.

$$R_{M+N}=\frac{D_M}{1-U_M{U}_N}{R}_N+R_M+\frac{D_M{U}_N{T}_M}{1-U_M{U}_N}$$

(19)

$U$ and $D$ are the backscattering and forward scattering ratios of the particle string, respectively. $R$ and $T$ are the reflectance and transmittance of the particle string, respectively. $R_M$ is equal to $T_M$ theoretically. We finally obtain:

$$R_{M+N}=\frac{D_M}{1-U_M{U}_N}{R}_N+\left(1+\frac{D_M{U}_N}{1-U_M{U}_N}\right)R_M$$

(20)

If the radiation source is switched to the right of the system, the subscript M and N in Equation (20) will be replaced with each other but $R_{M+N}$ remains the same.

$$R_{M+N}=\frac{D_N}{1-U_M{U}_N}{R}_M+\left(1+\frac{D_N{U}_M}{1-U_M{U}_N}\right)R_N$$

(21)

From the two equations above, we can obtain:

$$\frac{R_M}{R_N}=\frac{1-U_M{U}_N+D_N{U}_M-D_M}{1-U_M{U}_N+D_M{U}_N-D_N}$$

(22)

where equations of $U_N$ and $D_N$ can be found in Nomenclature. Since this is one-dimension case, we need to reassess the directional change coefficients for particle backscattering ratio $x_u$ and lateral scattering ratio $x_a$, they are denoted as $p_s$ and $p_{sa}$ (see Appendix A). If $M=1$, then $U_M=p_s{x}_u$, $R_M=\left[\left(1-p_a\right)x_a+\left(1-p\right)\left(x_u+x_d\right)\right]/2$, $D_M=p_s{x}_d$.Directional change coefficients for two-dimension case ($p$ and $p_a$) instead of $p_s$ and $p_{sa}$ are used to calculate $R_M$ because we need to exclude part of the laterally-scattered radiation. After taking $U_M$, $R_M$ and $D_M$ into Equation (22) we can get the form of $R_N$.

Now, let us consider the case in which the direction of incident radiation is perpendicular to the length of the particle string. Suppose there are $L$ particles in the string, light is incident on the $a$th left particle (Figure 8a). The substring left to this particle can be regarded as a block, the substring right to this particle can be regarded as another block (Figure 8b). When studying the radiation transfer within the particle string, only Components 2 and 3 need to be considered. Figure 8c,d give a simplified graph of radiation transfer within the particle string for these two components respectively. The reflectance for Component 2 and 3 can be derived like $R_{M+N}$, actually, we just need to replace $D$ with $U$ in Equation (20).

$$R_{cp2}=\frac{1}{1-U_l{U}_{r+1}}{R}_l+\frac{U_l}{1-U_l{U}_{r+1}}{R}_{r+1}$$

(23)

$$R_{cp3}=\frac{1}{1-U_r{U}_{l+1}}{R}_r+\frac{U_r}{1-U_r{U}_{l+1}}{R}_{l+1}$$

(24)

where $l$ and $r$ denote the numbers of particles to the left and right of the entrance particle, respectively and their values can be found in Figure 8a. Equations of $U_r$ and $R_r$ can be found in Nomenclature. Therefore, the final equation for global reflectance is:

$$R= m_e+\left(1-m_e\right)\tau \left[x_u+\frac{p_a{x}_a}{2}\left(R_{cp2}+R_{cp3}\right)\right]$$

(25)

The global transmittance is:

$$T= \left(1-m_e\right)\tau \left[x_d+\frac{p_a{x}_a}{2}\left(R_{cp2}+R_{cp3}\right)\right]$$

(26)

Six parameters are required to calculate the reflectance and transmittance of a particle string: total number of particles, location of the entrance particle, the refractive index, the diametrical absorbance, the particle backscattering and forward scattering ratios. If the location of the entrance particle moves from left to right in Figure 8a, we will first observe an increase of reflectance, then the reflectance reaches a plateau followed by a decline (Figure 9). The transmittance shows the same pattern. In this study, we select the plateau values as the final reflectance and transmittance.

4.4. Differences between LIBERTY and PROSPECT

The mechanisms of the forward simulation of LIBERTY and PROSPECT are shown in Figure 10. The comparisons between these two models are listed in

Table 3. Comparisons between Leaf Incorporating Biochemistry Exhibiting Reflectance and Transmittance Yields (LIBERTY) and The optical PROperties SPECTra model (PROSPECT5).

	PROSPECT5	LIBERTY	Note
Number of input parameters	6	9
Number of input biochemical parameters	5	6
Number of input structural parameters	1	3
Incorporated absorbance	kd	$\alpha d$	Equivalent
Sublayer morphology	Plate model	Compact particle model
Method to calculate global reflectance and transmittance	Stokes theory	Benford theory	Essentially the same
Boundary	Not considered	Not considered	same
Special treatment of incident radiation	Yes	No

. Despite the various chemical constituents which are measured in different bases, there is no substantial difference between these two models. Both models are built on a multiple plate model since their methods to calculate the global reflectance and transmittance (Benford theory and Stokes theory) are essentially the same. The main difference between LIBERTY and PROSPECT is the sublayer morphology. The sublayer is a flat plate in PROSPECT while in LIBERTY it is composed of compact particles. The distinct sublayer morphologies result in different methods to calculate the sublayer reflectance and transmittance. When calculating the sublayer transmittance, LIBERTY uses an indirect method. It first calculates the infinite reflectance which is then used to obtain the sublayer reflectance together with the sublayer transmittance.

For LIBERTY, the input biochemical parameters, the particle diameter $d$ and the calibrated absorption coefficients $k$ can be incorporated into one parameter: the diametrical absorbance $\alpha d$ (Figure 10a), which is equivalent to the parameter $kd$ (Figure 10b) in PROSPECT. However, the difference between $\alpha d$ and $kd$ is that $d$ can be eliminated in $kd$ because the contents of biochemical constituents can be measured in area basis for broadleaves. The detailed explanations are given as follows: According to Beer–Lambert law, the absorbance $kd$ would be dimensionless, which means the unit of $k$ would be cm⁻¹ if the unit of $d$ is cm. $k$ is the linear combination of the contents of the volumetric biochemical constituents $C_{vi}$ (g∙cm⁻³) and their corresponding absorption coefficients $k_i$ (cm⁻²·g).

$$k= {{\displaystyle \sum }}^{\text{}}{k}_i{C}_{vi}$$

(27)

However, the input biochemical constituents are area-based, i.e., rthe unit is g∙cm⁻². If the leaf thickness is D (unit cm), the volume-based biochemical constituents can be expressed as a function of area-based biochemical constituents ($C_i$) and D:

$$C_{vi}=C_i/D$$

(28)

After taking Equation (28) into Equation (27), we can get:

$$kd={{\displaystyle \sum }}^{\text{}}{k}_i{C}_i\frac{d}{D}=\frac{{{\displaystyle \sum }}^{\text{}}{k}_i{C}_i}{N}$$

(29)

where N is the number of sublayers. Actually, the particle diameter $d$ in $\alpha d$ can also be eliminated using the method above. But unlike broadleaves, the biochemical constituents of needle leaves cannot be expressed in area basis since needle leaves are not necessarily flat and the area of a needle leaf is difficult to measure. Area-based biochemical constituents of needle leaves only can be measured using projected area like Di Vittorio [10].

5. Materials and Methods

5.1. Ranges of Input Parameters

To conduct sensitivity analysis and inter-model comparisons, ranges of input parameters would be determined so as to represent more genuine physical and biochemical realities. There are five input parameters in this study: the relative refractive index $n$, the particle backscattering ratio $x_u$, the particle forward scattering ratio $x_d$, the diametrical absorbance $\alpha d$ and the structural parameter $N$. The biochemical constituents such as leaf chlorophyll content, leaf water content, lignin/cellulose content and nitrogen content have been incorporated into $\alpha d$; we did not consider the changes of specific biochemical constituents because: (1) the main limitations of LIBERTY exist in physical rather than biochemical aspects; (2) It is more easy to compare LIBERTY with other models such as PROSPECT which define different sets of biochemical components in various units.

The ranges of the refractive index n and diametrical absorbance $\alpha d$ were determined by the leaf optical properties model PROSPECT. PROSPECT was built by Jacquemoud and Baret [11], whereafter it underwent several improvements including recalibrations of refractive index and specific absorption coefficients [12]. In this study, the improved version PROSPECT5 was used (the source code can be downloaded from here). In order to get the ranges of the diametrical absorbance $\alpha d$, we first calculate the absorption coefficient of one elementary layer $k\left(\lambda \right)$ which is a linear combination of the content of leaf constituents $C_i$ and their corresponding specific absorption coefficients $k_i\left(\lambda \right)$:

$$k\left(\lambda \right)={\displaystyle \sum }_i{k}_i\left(\lambda \right)\frac{C_i}{N}$$

(30)

where $N$ is the structural parameter, and $k_i\left(\lambda \right)$ was taken from PROSPECT5 while $C_i$ and $N$ were read from the LOPEX93 dataset, which was downloaded from here. A description of this dataset can be found here. The range of the diametrical absorbance $\alpha d$ is $\left[\min \left(k\right), \max \left(k\right)\right]$. The range of structural parameter $N$ can be read from this dataset as well.

The determination of the $x_u$ range is more difficult. According to the study of Mandelis, Boroumand and Bergh [7], the value of $x_u$ would lie in [0.1, 0.4]. Simonot, et al. [13] also calculated the particle backscattering ratio using ray tracing technique and found the value is in the range 0.25–0.35. However, in their study, all upper hemispherical fluxes were counted as backscattered radiation which thus included part of the laterally-scattered radiation. In this study, we assume $x_u$ is located in [0, 0.5]. A statistic summary of parameter ranges is listed in Table 1.

5.2. Sensitivity Analysis

Sensitivity analysis (SA) helps to figure out contributions of variability in input parameters to the global variability in model outputs. It has been widely used in mathematical models with complicated input parameters. Generally, there are two types of SA: local SA and global SA. Their major difference is whether consider the interactions between model parameters. In global SA, both individual parameter and interactions between parameters are supposed to explain the global variability while in local SA, only an individual parameter is considered. The theory of global SA can be found in previous literature such as by Ceccato, et al. [14] and Wang, et al. [15]. Suppose a model needs three input parameters. The total variance can be explained by the variance of single parameter and the variance of interactions between parameters:

$$V=V_1+V_2+V_3+V_{12}+V_{13}+V_{23}+V_{123}$$

(31)

where $V_i$ denotes the variance of parameter $i$, $V_{ij}$ is the variance of the interaction between parameters $i$ and $j$, and $V_{ijk}$ is the variance of the interaction between parameters $i$, $j$ and $k$.

The first-order sensitivity index $S_i$ for parameter $i$ is defined as the ratio between partial input variance of $i$ and the total variance. It represents the percentage of the output variance that is accounted for by input parameter $i$.

$$S_i= V_i/V$$

(32)

The total-order sensitivity index $S{T}_i$ for parameter $i$ is defined as the ratio of total partial input variance of $i$ (including single parameter and interactions between parameter $i$ and other parameters) to the total variance.

$$S{T}_i=(V_i+{\displaystyle \sum }_j{V}_{ij}+V_{ijk})/V$$

(33)

The difference between $S{T}_i$ and $S_i$ is the percentage of the output variance that is accounted for by interactive effects of parameter $i$ with others.

A Matlab software tool—global sensitivity analysis toolbox (GSAT) [16] was used in this study to perform global sensitivity analysis. The source code can be found here. The original program can only produce samples when input parameters are independent to each other. However, in LIBERTY $x_u$ and $x_d$ are constrained under $1\le x_u+x_d\le 1$. In order to meet this special need, we made some modifications over the source code [17,18]. The specific method is to define two new variables $A$ and $B$ whose values are:

$$A=x_u+x_d$$

(34)

$$\mathrm{B}=x_u/\left(x_u+x_d\right)$$

(35)

Both $A$ and $B$ lie in [0, 1]. $A$ and $B$ are regarded as input parameters instead of $x_u$ and $x_d$. After samples have been produced, the values of $x_u$ and $x_d$ can be derived by the two equations above. There are two sensitivity analysis methods commonly used in previous literature: the Fourier amplitude sensitivity test (FAST) [19] and the Sobel method [20,21]. The Sobel method was adopted for inter-model comparisons while FAST was used for the sensitivity analysis of the particle string model. A quasi-monte Charlo algorithm is employed by GSAT, which produces two sets of samples: one set for evaluation and the other for replacements of complementary variables. The constraint of $1\le x_u+x_d\le 1$ cannot be ensured during the replacements.

5.3. Inter-Model Comparisons

The original LIBERTY, PROSPECT5 and our improved version of LIBERTY (LIBERTY_im) were compared with each other in this study. Despite the variety of input biochemical constituents, the main difference between LIBERTY and PROSPECT is in their sublayer morphologies (see Section 4.4). The absorbance kd used in PROSPECT is equivalent to the diametrical absorbance $\alpha d$ in LIBERTY. The comparisons were conducted in three aspects: the sublayer backscattering ratio, the sublayer reflectance and transmittance, the global reflectance and transmittance. Since the numbers of the input parameters for three models are inconsistent (

Table 4. Inputs for three models *.

	Models
	LIBERTY	LIBERTYim **	PROSPECT5
Outputs	Inputs
$X_u{X}_d{X}_u$ and $X_d$	$n,\alpha d,x_u$	$n,\alpha d,x_u$	$n,\alpha d$
$R_s{T}_s{R}_s$ and $T_s$	$n,\alpha d,x_u$	$n,\alpha d,x_u$	$n,\alpha d$
$RTR$ and $T$	$n,\alpha d,x_u,N$	$n,\alpha d,x_u,N$	$n,\alpha d,N$

* the definitions of the parameters are given in Nomenclature. ** improved LIBERTY in this study.

), the numbers of sensitivity indexes are inconsistent as well. There is no scattering ratio in PROSPECT, so the comparisons of sublayer scattering ratios were conducted only between LIBERTY and LIBERTY_im. The particle backscattering ratio $x_u$ was assumed to be equal to the forward scattering ratio $x_d$ in both LIBERTY and LIBERTY_im for better comparisons between the two models. As a result, the sublayer backscattering and forward scattering ratios were equal to each other as well. The refractive index is fixed in PROSPECT5, but in this study, it serves as a variable.

6. Results

6.1. Sublayer Scattering Ratios Considering the Directional Changes of the Particle Scattering Ratios

In general, the sensitivity patterns of the sublayer backscattering ratio $X_u$ for LIBERTY and LIBERTY_im were similar (Figure 11a). Among three input parameters, the particle backscattering ratio $x_u$ explained the most of the variances of $X_u$ for both models. The second sensitive parameter was the diametrical absorbance $\alpha d$. Neither LIBERTY nor LIBERTY_im was sensitive to the refractive index $n$ and the effects of interactions among parameters were insignificant. Although similar in sensitivity patterns, the sublayer backscattering ratios of LIBERTY and LIBERTY_im were quite different in values if the same input parameters were entered (Figure 11b).

6.2. The Sublayer Reflectance and Transmittance

The important factor for the sublayer reflectance and transmittance of all three models was the diametrical absorbance $\alpha d$ (Figure 12). Neither LIBERTY nor LIBERTY_im was sensitive to the refractive index $n$.

The sublayer reflectance $R_s$ of LIBERTY_im was more sensitive to the diametrical absorbance $\alpha d$ but less sensitive to the particle backscattering ratios $x_u$ than LIBERTY. In contrast, the global sensitivity index of the sublayer transmittance $T_s$ to $\alpha d$ decreased slightly after including the directional changes of the particle scattering ratios and the percentage of output variance explained by interactions among parameters changed negligibly. Sensitivity of the sublayer transmittance to the particle backscattering ratios $x_u$ did not change too much. The sublayer reflectance and transmittance of LIBERTY_im had similar patterns of sensitivities to input parameters. Among the three models, PROSPECT needs the least number of parameters to calculate the sublayer reflectance and transmittance (the refractive index n and the diametrical absorbance $\alpha d$ only). Interactions between these two parameters could explain a great deal of output variance of the sublayer reflectance of PROSPECT. Compared with the sublayer transmittance, the sublayer reflectance of PROSPECT is much more sensitive to the refractive index $n$. The calculated sublayer reflectance and transmittance of LIBERTY and LIBERTY_im were quite different if the same input parameters were entered (Figure 13a).

6.3. The Global Reflectance and Transmittance

Generally, the diametrical absorbance $\alpha d$ was the dominant factor contributing to both global reflectance and transmittance of all three models while the refractive index $n$ was the insignificant factor (Figure 14). The global reflectance of LIBERTY was sensitive to the particle backscattering ratio $x_u$ while the interactions among input parameters was not very important. Meanwhile, the number of sublayers $N$ and the interactions among input parameters played important roles in the output variance of the global transmittance of LIBERTY. For LIBERTY_im, the global reflectance was insensitive to other parameters apart from the diametrical absorbance $\alpha d$. However, compared with the global reflectance, the global transmittance was much more sensitive to the number of sublayers $N$ and the interactions among input parameters. For PROSPECT5, parameters except the diametrical absorbance $\alpha d$ had little impact on the variance of the global reflectance.

The sensitivity pattern of the global reflectance and transmittance of LIBERTY was quite different from that of PROSPECT. But after being physically corrected, LIBERTY had similar pattern to PROSPECT in simulated global reflectance and transmittance. Considering the directional changes of the scattering ratios witnessed an increased sensitivity of the global reflectance and transmittance to the diametrical absorbance. However, the sensitivities of the global reflectance to the particle backscattering ratio $x_u$, the number of sublayers $N$ and the global transmittance to interactions among parameters decreased sharply. When the same input parameters were entered into LIBERTY and LIBERTY_im, great differences in the infinite reflectance, the sublayer and global reflectance and transmittance were observed (Figure 13).

6.4. Reflectance and Transmittance of a Particle String

The differences of sensitivity between the reflectance and transmittance are quite small (see Figure 15). Both of them are most sensitive to the diametrical absorbance $\alpha d$ followed by the interactions among parameters. Other factors such as the number of particles $N$, the refractive index n, the particle scattering ratios $x_u$ and $x_d$ contribute little to output reflectance and transmittance.

7. Discussion

The directional changes of the particle scattering ratios have great impact on the sublayer backscattering ratios (Figure 11b). The sublayer reflectance and transmittance were calculated based on the sublayer scattering ratios with consideration of the entrance particle. Since the sublayer backscattering ratio and forward scattering ratio were assumed to be equal in this study, the sensitivity differences of the sublayer reflectance and transmittance to input parameters between LIBERTY and LIBERTY_im are due to the distinct treatments of the radiation components related to the entrance particle. As explained in Section 4.4, the methods used by LIBERTY and PROSPECT to calculate the global reflectance and transmittance are essentially the same. However, the sensitivity analysis of the global reflectance and transmittance showed different patterns for these two models. By contrast, LIBERTY_im has similar sensitivity patterns to PROSPECT with respect to the global reflectance and transmittance. The results demonstrate indirectly the limitations of LIBERTY and the validity of LIBERTY_im. The advantage of LIBERTY_im is that neither its global reflectance nor transmittance is sensitive to the particle scattering ratios.

Compared with PROSPECT, LIBERTY needs an additional parameter (i.e., the particle backscattering ratio $x_u$) to obtain the sublayer reflectance and transmittance, which is determined by its distinct sublayer morphology. The particle forward scattering ratio $x_d$ was assumed to be equal to the particle backscattering ratio $x_u$. However, the condition for this equality (i.e., $\alpha d<1$) does not always hold (see Table 1), indicating independence of $x_d$. Thus, $x_d$ would be entered into LIBERTY as an input parameter under the constraint of $1\le x_u+x_d\le 1$. The introductions of the additional parameters will result in more uncertainties. Actually, the sublayer scattering ratios and the diametrical absorbance are related with each other in an unknown fashion. Once their relationship can be determined by mathematical equations, LIBERTY will have fewer parameters. The main limitation of LIBERTY is the assumption that the sublayer is infinite extending. Needle leaves are narrow and thin, hence, they cannot be regarded as endless in width. If the boundary constraints of needle leaves are not considered, LIBERTY is more suitable for modelling the optical properties of broadleaves instead of needle leaves.

In this study, we worked out the reflectance and transmittance of a particle string, which is a theoretical model. The particle string model may be capable of modelling leaf optical properties for needle leaves. If the model is used for this purpose, the particle in this model will no longer denote the cell anymore; the whole needle is regarded as a particle string. The greatest advantage of this particle string model is the inclusion of the special morphology of needle leaves. The equivalent diameter $d$ of particle can be calculated using measured thickness and width. The particle number can be assumed to be infinite since the needle leaf is long enough, thus the number of particles N can be eliminated. There are only four parameters left: refractive index $n$, diametrical absorbance $\alpha d$, particle backscattering and forward scattering ratios ($x_u$ and $x_d$). $x_u$ and $x_d$ can be treated as parameters controlling the anatomical morphology of needle leaves. However, there are several issues need to be addressed if the particle string model can be used. Firstly, the absorption coefficients $\alpha$ may not reflect chemical reality, it is just a calibrated parameter used in the model, but this does not matter since the ultimate goals are to simulate the reflectance and transmittance and to estimate biochemical concentrations. The intermediate products such as absorption coefficients will not affect the performance of the model. Secondly, $x_u$ and $x_d$ are related to $\alpha d$ and $n$. If we regard the particle scattering ratios and the diametrical absorbance as independent variables, the physical reality of the particle scattering is violated. However, it seems hard to determine their relationships. Once their relations can be expressed precisely by mathematical equations, the parameters $x_u$ and $x_d$ can be eliminated. There are two parameters left: refractive index $n$ and diametrical absorbance $\alpha d$. Anatomical morphology of needle leaves can be controlled by diametrical absorbance $\alpha d$. Thus the parsimony of this model needs to be checked. Lastly, the particle ignores the intercellular air spaces inside the needle leaves. LIBERTY describes the intercellular organizations inside the leaf in a more delicate manner than PROSPECT. However, it cannot be used to model optical properties for needle leaves. Both LIBERTY and PROSPECT has considered intercellular air spaces. To include intercellular air spaces into the particle string model, differences between broadleaves and needle leaves would be noticed. Although this model is theoretically capable of modeling optical properties for needle leaves, experiments are needed to test its applicability.

8. Conclusions

Our study pointed out several limitations existing in LIBERTY according to clues of these limitations. Based on the findings, theoretical improvements of LIBERTY were given. The improved version of LIBERTY (LIBERTY_im) has limited differences to PROSPECT except in their distinct sublayer morphologies. Unless the special realities of needle leaves (boundary constraints) are included, LIBERTY cannot be truly used for modelling optical properties for needle leaves. However, it is difficult to introduce the boundary constraints into a sublayer in which a particle denotes a cell. However, the reflectance and transmittance of a string composed limited number of particles can be calculated. We worked out the reflectance and transmittance of a particle string which is composed of a limited number of particles. The particle string model has great potential for needle leaf optical properties modelling. It includes leaf morphology and radiation scattered in all directions around the needle leaves can be calculated using this particle string. However, delicate experiments are needed to test the theory.