Development of Analytical Model to Describe Reflectance Spectra in Leaves with Palisade and Spongy Mesophyll

Abstract

Remote sensing plays an important role in plant cultivation and ecological monitoring. This sensing is often based on measuring spectra of leaf reflectance, which are dependent on morphological, biochemical, and physiological characteristics of plants. However, interpretation of the reflectance spectra requires the development of new tools to analyze relations between plant characteristics and leaf reflectance. The current study was devoted to the development, parameterization, and verification of the analytical model to describe reflectance spectra of the dicot plant leaf with palisade and spongy mesophyll layers (on the example of pea leaves). Four variables (intensities of forward and backward collimated light and intensities of forward and backward scattered light) were considered. Light reflectance and transmittance on borders of lamina (Snell’s and Fresnel’s laws), light transmittance in the palisade mesophyll (Beer–Bouguer–Lambert law), and light transmittance and scattering in the spongy mesophyll (Kubelka–Munk theory) were described. The developed model was parameterized based on experimental results (reflectance spectra, contents of chlorophylls and carotenoid, and thicknesses of palisade and spongy mesophyll in pea leaves) and the literature data (final R² was 0.989 for experimental and model-based reflectance spectra). Further model-based and experimental investigations showed that decreasing palisade and spongy mesophyll thicknesses in pea leaves (from 35.5 to 25.2 µm and from 58.6 to 47.8 µm, respectively) increased reflectance of green light and decreased reflectance of near-infrared light. Similarity between model-based and experimental results verified the developed model. Thus, the model can be used to analyze leaf reflectance spectra and, thereby, to increase efficiency of the plant remote and proximal sensing.

1. Introduction

Terrestrial plants can be affected by numerous environmental factors influencing their growth and development, productivity, physiological processes, and content of biochemical compounds. Early-revealing changes in the characteristics of plants play an important role in their agricultural cultivation and ecological monitoring and require the development of methods of remote sensing [1,2]. Optical methods, including RGB imaging, fluorescence imaging, thermal imaging, and imaging of spectral characteristics of reflected light, are effective tools for plant remote sensing [3,4] because they are non-invasive, relatively simple, and fast.

Particularly, multispectral imaging, which is based on measuring reflectance in several narrow spectral bands, and hyperspectral imaging, which is based on measuring reflectance spectra, are widely used for plant remote sensing [3,4]. It is known that leaf reflectance in narrow spectral bands can be related to specific plant characteristics, including the content of photosynthetic pigments [5,6,7], water content [8,9,10], nitrogen content [11,12], leaf area index and biomass [12,13,14,15], photosynthetic processes and their stress changes [16,17,18,19,20], etc. Calculation of reflectance indices, which are mainly dimensionless indicators based on reflectance at two or three wavelengths, is used to improve these relationships and to increase efficiency of estimation of plant characteristics [3,4].

However, relationships of specific reflectance indices in leaves to specific plant characteristics can be varied [3]; for example, this variability was shown [21,22,23] for photochemical reflectance index based on reflectance at 531 and 570 nm. It is known that this index is dependent on the acidification of the chloroplast lumen (through transitions in the xanthophyll cycle [17,19] and through the chloroplast shrinkage [20]) and, therefore, should be potentially related to photosynthetic processes. Moreover, our earlier meta-analysis of experimental studies [23] showed that average linear correlation coefficients between values of photochemical reflectance indices and photosynthetic parameters (the photosynthetic light use efficiency, potential quantum yield of photosystem II, and non-photochemical quenching of chlorophyll fluorescence) are about 0.5–0.6.

There are several possible reasons for the variability of relationships of reflectance at specific wavelengths (and, therefore, reflectance indices) to plant characteristics [3]. First, a photosynthetic pigment composition, which can be dependent on plant species [24,25], individual development [11,26,27,28], and environmental conditions [24], strongly influences reflectance spectra (through changes in spectra of the light absorption in the leaf lamina). Second, differences in the leaf anatomy can also increase variability of reflectance spectra [3]. It is known that reflectance of visible light and near-infrared light (NIR) is different on adaxial and abaxial leaf surfaces [29,30,31], which is caused by different optical properties of palisade and spongy mesophylls in dicot plants [29,32]. Reflectance is also dependent on the leaf thickness [25,32,33] and presence of hairs [34]. Third, leaf orientation (angle between the leaf surface and incident light) [25,29] and their fluttering under the wind [35] can also modify reflectance spectra. Fourth, changes in phenological stages of plants and their senescence strongly influence reflectance spectra [6,26,36,37]. The influence can have complex reasons, including modification of content of photosynthetic pigments and leaf orientation [37], that is, it can be caused by mechanisms, which are noted above.

Thus, variability of leaf reflectance spectra, which decreases relationships of parameters of these spectra to plant characteristics, can limit analysis of results of the plant remote sensing and requires the development of tools to minimize these limitations. Mathematical models of optical characteristics of leaves, which consider content of photosynthetic pigments and anatomical structure, can be potentially used as these tools. In accordance with [38], there are several types of models to describe the optical properties of leaves, namely, Beer–Bouguer–Lambert law-based models, plate models, N-flux models (or Kubelka–Munk theory-based models), compact spherical particle models, radiative transfer theory-based models, stochastic models, and ray tracing models.

Kubelka–Munk theory-based models, which are widely used to analyze light reflectance and transmittance in leaves and canopy [39,40,41,42,43,44,45,46], are relatively simple (analytical solutions can be derived) and often provide accurate descriptions of the optical properties of plants. The models consider two (forward and backward) [39,41,43,44,45,46] or four (forward collimated, forward scattered, backward collimated, and backward scattered) [40,42] light flows. Particularly, the Kubelka–Munk theory-based model can be used to derivate light absorption coefficients from reflectance spectra of canopy [44] and leaf [45,46,47] or to analyze light-dependent chloroplast movements [48]. It should also be noted that description of plant leaf can be difficult due to several layers with different optical properties (mainly, palisade and spongy mesophyll layers, which are typical for dicot plants) [38]. There is the Kubelka–Munk theory-based model (with forward and backward light flows) [41] describing different layers in plant leaves; this study shows that Kubelka–Munk theory-based models can be potentially used to describe reflectance of leaf with different layers.

Moreover, simple Beer–Bouguer–Lambert law-based models can be used at weak light scattering [38]. Considering the low light scattering coefficient in the leaf palisade mesophyll [49], these models are probable to be effective to describe the optical properties of the layer. Using Beer–Bouguer–Lambert law-based models for the palisade mesophyll should simplify an analytical description of light reflectance in leaves.

Thus, the current study was devoted to the development, parameterization, and verification of an analytical model to describe leaf reflectance spectra in dicot higher plants (with palisade and spongy mesophyll layers), which was based on both the Kubelka–Munk theory and the Beer–Bouguer–Lambert law as well as on other widely used approaches (particularly the Snell’s and Fresnel’s laws). The model can be potentially used as the tool for analysis of results of the plant remote sensing on the basis of reflectance measurements and for the development of new indicators to estimate plant characteristics.

2. Description of the Mathematical Model of Light Reflectance in Leaves

2.1. Basic Structure and Variables of the Model

In accordance with Figure 1, the leaf was assumed to be an optical system, including two main optical layers, namely, the palisade mesophyll layer with high light absorption ($a_P$) and low light scattering ($s_P$) coefficients and the spongy mesophyll layer with high light absorption ($a_{Sp}$) and high light scattering ($s_{Sp}$) coefficients [49]. Epidermal layers had small thickness and weak color (Figure 2a) and were not considered in the model. Borders “air-leaf” and “leaf-air” were additionally described in the model. Four variables were considered in the current study (Figure 1), including intensities of the forward collimated light (I_C), forward scattered light (I_S), backward collimated light (J_C), and backward scattered light (J_S) [40,42].

2.2. Equations Describing Light Reflectance and Transmittance on Borders “Air-Leaf” and “Leaf-Air”

It was assumed that there were two types of transmittances of the collimated light across each border (“air-leaf” and “leaf-air”): (i) across a smooth surface and (ii) across a rough surface. The fraction of the rough surface (F_S) was the model parameter; the fraction of the smooth surface was calculated as 1-F_S.

Relationships between the angles of incidence and refraction were described by Equations (1) and (2), which were based on Snell’s law:

$${\beta }_{I1}=\arcsin \left({\displaystyle \frac{n_O}{n_I}}\sin \left({\beta }_{O1}\right)\right)$$

(1)

$${\beta }_{I2}=\arcsin \left({\displaystyle \frac{n_O}{n_I}}\sin \left({\beta }_{O2}\right)\right)$$

(2)

where n_O and n_I are the refractive indices in air (n_O = 1) and in leaf (n_I = 1.415 [50]); β_O1 and β_O2 are the angles of incidence on adaxial and abaxial leaf surfaces, respectively; and β_I1 and β_I2 are the angles of refraction under adaxial and abaxial leaf surfaces, respectively.

Fresnel’s law-based Equations (3) and (4) were used for the calculation of transmittance coefficients for the collimated light transfer from air to leaf across the smooth surfaces [51,52]:

$${T_{Ic}}^{OI}=1-{\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{\sin \left({\beta }_{O1}-{\beta }_{I1}\right)}{\sin \left({\beta }_{O1}+{\beta }_{I1}\right)}}\right)}^2+{\left({\displaystyle \frac{\tan \left({\beta }_{O1}-{\beta }_{I1}\right)}{\tan \left({\beta }_{O1}+{\beta }_{I1}\right)}}\right)}^2\right]$$

(3)

$${T_{Jc}}^{OI}=1-{\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{\sin \left({\beta }_{O2}-{\beta }_{I2}\right)}{\sin \left({\beta }_{O2}+{\beta }_{I2}\right)}}\right)}^2+{\left({\displaystyle \frac{\tan \left({\beta }_{O2}-{\beta }_{I2}\right)}{\tan \left({\beta }_{O2}+{\beta }_{I2}\right)}}\right)}^2\right]$$

(4)

where T_Ic^OI and T_Jc^OI are transmittance coefficients for the collimated light transfer from air to leaf on adaxial and abaxial surfaces, respectively. Thus, Equations (5) and (6) were used for the calculation of the intensity of the collimated light transmittance into the leaf across adaxial and abaxial surfaces (I_C(0) and J_C(h + l), respectively):

$$I_C(0)=I_0\cdot \left(1-F_S\right)\cdot {T_{Ic}}^{OI}$$

(5)

$$J_C(h+l)=J_0\cdot \left(1-F_S\right)\cdot {T_{Jc}}^{OI}$$

(6)

where I₀ and J₀ are the intensities of the forward and backward collimated light directed to adaxial and abaxial leaf surfaces (the incident light), respectively; h is the thickness of the palisade mesophyll layer; and l is the thickness of the spongy mesophyll layer.

Equations (7) and (8) were used for the calculation of the intensity of the collimated light reflecting from adaxial and abaxial leaf surfaces in air (J_C^RO and I_C^RO, respectively):

$${J_C}^{RO}=I_0\cdot \left(1-F_S\right)\cdot \left(1-{T_{Ic}}^{OI}\right)$$

(7)

$${I_C}^{RO}=J_0\cdot \left(1-F_S\right)\cdot \left(1-{T_{Jc}}^{OI}\right)$$

(8)

Equations (9) and (10) were used for the calculation of the transmittance coefficients for collimated light transfer from leaf to air:

$${T_{Ic}}^{IO}=1-{\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{\sin \left({\beta }_{I1}-{\beta }_{O1}\right)}{\sin \left({\beta }_{I1}+{\beta }_{O1}\right)}}\right)}^2+{\left({\displaystyle \frac{\tan \left({\beta }_{I1}-{\beta }_{O1}\right)}{\tan \left({\beta }_{I1}+{\beta }_{O1}\right)}}\right)}^2\right]$$

(9)

$${T_{Jc}}^{IO}=1-{\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{\sin \left({\beta }_{I2}-{\beta }_{O2}\right)}{\sin \left({\beta }_{I2}+{\beta }_{O2}\right)}}\right)}^2+{\left({\displaystyle \frac{\tan \left({\beta }_{I2}-{\beta }_{O2}\right)}{\tan \left({\beta }_{I2}+{\beta }_{O2}\right)}}\right)}^2\right]$$

(10)

where T_Ic^IO and T_Jc^IO are the transmittance coefficients for the collimated light transfer from leaf to air on abaxial and adaxial leaf surfaces, respectively. Equations (9) and (10) could not be solved at β_I1 and β_I2, which were more than about 45° because light was fully reflected in this case; the case was considered to describe transmittance and reflectance of the scattered light. It should be additionally noted that T_Jc^IO = T_Jc^OI and T_Ic^IO = T_Ic^OI for the collimated light.

Equations (11) and (12) were used for the calculation of the intensity of the collimated light transferring from leaf to air across adaxial and abaxial surfaces (J_C^T and I_C^T, respectively):

$${J_C}^T=J_C\left(0\right)\cdot \left(1-F_S\right)\cdot {T_{Jc}}^{IO}$$

(11)

$${I_C}^T=I_C\left(h+l\right)\cdot \left(1-F_S\right)\cdot {T_{Ic}}^{IO}$$

(12)

Equations (13) and (14) were used for the calculation of the intensity of the collimated light reflecting from adaxial and abaxial leaf surfaces in the lamina (I_C^RI and J_C^RI, respectively):

$${I_C}^{RI}=J_C\left(0\right)\cdot \left(1-F_S\right)\cdot \left(1-{T_{Jc}}^{IO}\right)$$

(13)

$${J_C}^{RI}=I_C\left(h+l\right)\cdot \left(1-F_S\right)\cdot \left(1-{T_{Ic}}^{IO}\right)$$

(14)

Transmittance coefficients for the collimated light transfer across the rough surfaces and for the scattered light transfer across both smooth and rough surfaces were described by Equations (15)–(18):

$${T_{Is}}^{OI}={\displaystyle \frac{2}{\pi }}{\displaystyle \underset{0}{\overset{\pi /2}{\int }}{T_{Ic}}^{OI}}d{\beta }_{O1}$$

(15)

$${T_{Js}}^{OI}={\displaystyle \frac{2}{\pi }}{\displaystyle \underset{0}{\overset{\pi /2}{\int }}{T_{Jc}}^{OI}}d{\beta }_{O2}$$

(16)

$${T_{Is}}^{IO}={\displaystyle \frac{2}{\pi }}{\displaystyle \underset{0}{\overset{\pi /2}{\int }}{T_{Ic}}^{IO}}d{\beta }_{I1}$$

(17)

$${T_{Js}}^{IO}={\displaystyle \frac{2}{\pi }}{\displaystyle \underset{0}{\overset{\pi /2}{\int }}{T_{Jc}}^{IO}}d{\beta }_{I2}$$

(18)

where T_Is^OI and T_Js^OI are the transmittance coefficients for the light transfer from air to leaf on adaxial and abaxial surfaces, respectively, and T_Is^IO and T_Js^IO are transmittance coefficients for the light transfer from leaf to air on abaxial and adaxial surfaces, respectively. It should be noted that T_Is^OI = T_Js^OI and T_Is^IO = T_Js^IO; thus, T_s^OI, equaling to T_Is^OI (=T_Js^OI) and T_s^IO equaling to T_Is^IO (=T_Js^IO) were used in the analysis. On basis of numerical calculation, it was shown that T_s^OI ≈ 0.866 and T_s^IO ≈ 0.469. These values were used for modeling.

Equations (19) and (20) were used for the calculation of the intensity of the scattered light transferring from air to leaf across adaxial and abaxial surfaces (I_S(0) and J_S(h + l), respectively):

$$I_S(0)=I_0\cdot F_S\cdot {T_s}^{OI}$$

(19)

$$J_S(h+l)=J_0\cdot F_S\cdot {T_s}^{OI}$$

(20)

Equations (21) and (22) were used for the calculation of the intensity of the scattered light reflecting from adaxial and abaxial leaf surfaces in air (J_S^RO and I_S^RO, respectively):

$${J_S}^{RO}=I_0\cdot F_S\cdot \left(1-{T_s}^{OI}\right)$$

(21)

$${I_S}^{RO}=J_0\cdot F_S\cdot \left(1-{T_s}^{OI}\right)$$

(22)

Equations (23) and (24) were used for the calculation of the intensity of the scattered light transferring from leaf to air across adaxial and abaxial surfaces (J_S^T and I_S^T, respectively):

$${J_S}^T=J_C\left(0\right)\cdot F_S\cdot {T_s}^{IO}+J_S\left(0\right)\cdot {T_s}^{IO}$$

(23)

$${I_S}^T=I_C\left(h+l\right)\cdot F_S\cdot {T_s}^{IO}+I_S\left(h+l\right)\cdot {T_s}^{IO}$$

(24)

Equations (25) and (26) were used for the calculation of the intensity of the scattered light reflecting from adaxial and abaxial leaf surfaces in the lamina (I_S^RI and J_S^RI, respectively):

$${I_S}^{RI}=J_C\left(0\right)\cdot F_S\cdot \left(1-{T_s}^{IO}\right)+J_S\left(0\right)\cdot \left(1-{T_s}^{IO}\right)$$

(25)

$${J_S}^{RI}=I_C\left(h+l\right)\cdot F_S\cdot \left(1-{T_s}^{IO}\right)+I_S\left(h+l\right)\cdot \left(1-{T_s}^{IO}\right)$$

(26)

Equations (27) and (28) were used for the calculation of the total intensities of light flows from adaxial and abaxial leaf surfaces to air (J_out¹ and I_out¹, respectively):

$${J_{out}}^1={J_C}^{RO}+{J_C}^T+{J_S}^{RO}+{J_S}^T$$

(27)

$${I_{out}}^1={I_C}^{RO}+{I_C}^T+{I_S}^{RO}+{I_S}^T$$

(28)

where “1” shows that these light intensities were calculated on the basis of the first iteration of the light propagation through leaf (see Section 2.5 for details).

2.3. Equations Describing Light Transmittance in the Palisade Mesophyll Layer

It is known that the palisade mesophyll layer has the low light scattering coefficient ($s_P$ = 5 cm⁻¹) in comparison to the spongy mesophyll layer ($s_{Sp}$ = 1000 cm⁻¹) [49]. It is known that the thickness of the palisade mesophyll layer (h) is often less than 70–80 μm [49,53,54,55]. It means that the probability of the primary light scattering in the palisade mesophyll layer is less than 3–4%; that is, a description of this scattering is not necessary to describe the light propagation through the layer.

As a result, we used Beer–Bouguer–Lambert law [38] as the basis of the description of the light propagation through the palisade mesophyll layer. Equations (29)–(32) were used for this description:

$$I_C(x)=I_C(0)\cdot e^{-{\displaystyle \frac{a_P}{\cos \left({\beta }_{I1}\right)}}x}$$

(29)

$$J_C(x)=J_C\left(h\right)\cdot e^{-{\displaystyle \frac{a_P}{\cos \left({\beta }_{I2}\right)}}\left(h-x\right)}$$

(30)

$$I_S(x)=I_S(0)\cdot e^{-2a_{P}x}$$

(31)

$$J_S(x)=J_S\left(h\right)\cdot e^{-2a_P\left(h-x\right)}$$

(32)

where x is the coordinate, $a_P$ is the light absorption coefficient, and “2” is the coefficient showing increase in the light absorption (and scattering) for the scattered light (on the basis of pathlength averaging over a hemisphere [42]).

Based on these equations, Equations (33)–(36), which described the light intensity on borders of the palisade mesophyll layer, were derived:

$$I_C(h)=I_C(0)e^{-{\displaystyle \frac{a_P}{\cos \left({\beta }_{I1}\right)}}h}$$

(33)

$$J_C(0)=J_C\left(h\right)\cdot e^{-{\displaystyle \frac{a_P}{\cos \left({\beta }_{I2}\right)}}h}$$

(34)

$$I_S(h)=I_S(0)\cdot e^{-2a_{P}h}$$

(35)

$$J_S(0)=J_S\left(h\right)\cdot e^{-2a_{P}h}+{J_S}^{Add}$$

(36)

where J_C(0) and I_C(h) are the intensities of collimated backward and forward light on upper and lower borders of the palisade mesophyll layer, respectively; J_S(0) and I_S(h) are intensities of scattered backward and forward light on upper and lower borders of the palisade mesophyll layer, respectively; and J_S^Add is the additional scattered light. J_S^Add is caused by the scattering (and changing light direction) of the scattered (I_S(x)) and collimated (I_C(x)) forward light in the palisade mesophyll layer. The J_S^Add should have low intensity; however, it can be important for the plant remote sensing based on measuring reflectance at red and light spectral bands, which have high light absorption.

We used Equation (37) to calculate J_S^Add:

$${J_S}^{Add}=I_C(0)\cdot {\displaystyle \frac{s_P\left(1-f\right)}{\cos \left({\beta }_{I1}\right)}}{\displaystyle \underset{0}{\overset{h}{\int }}{e}^{-\frac{a_P}{\cos \left({\beta }_{I1}\right)}x}\cdot e^{-2a_{P}x}}dx+2I_S(0)\cdot s_P\left(1-f\right){\displaystyle \underset{0}{\overset{h}{\int }}{e}^{-2a_{P}x}\cdot e^{-2a_{P}x}}dx$$

(37)

where $s_P$ is the light scattering coefficient in the palisade mesophyll layer and f is the asymmetry factor, which can describe the anisotropy of scattering (we assumed that f = 0.5; i.e., the asymmetry was absent). Equation (38) is the solution of Equation (37):

$${J_S}^{Add}=I_C(0)\cdot {\displaystyle \frac{s_P\left(1-f\right)}{a_P\left(1+2\cos \left({\beta }_{I1}\right)\right)}}\cdot \left(1-e^{-\left({\displaystyle \frac{a_P}{\cos \left({\beta }_{I1}\right)}}+2a_P\right)h}\right)+I_S(0)\cdot {\displaystyle \frac{s_P\left(1-f\right)}{2a_P}}\cdot \left(1-e^{-4a_{P}h}\right)$$

(38)

2.4. Equations Describing Light Transmittance and Scattering in the Spongy Mesophyll Layer

Considering the high light scattering coefficient in the spongy mesophyll layer ($s_{Sp}$ = 1000 cm⁻¹ [49]), we used the Kubelka–Munk model with four light flows [40,42] to describe optical properties in this layer. We used modified coordinate x₁ (x₁ = x − h) to simplify the analysis. x₁ can be used from h to h + l only, where h and l are the thicknesses of the palisade and spongy mesophyll layers, respectively. System of Equation (39) shows the initial Kubelka–Munk equations:

$$\begin{gathered} {\displaystyle \frac{d{I}_C(x_1)}{d{x}_1}}=-{\displaystyle \frac{a_{Sp}+s_{Sp}}{\cos \left({\beta }_{I1}\right)}}\cdot I_C(x_1) \\ \frac{d{J}_C(x_1)}{d{x}_1}={\displaystyle \frac{a_{Sp}+s_{Sp}}{\cos \left({\beta }_{I2}\right)}}\cdot J_C(x_1) \\ \begin{array}{l}{\displaystyle \frac{d{I}_S(x_1)}{d{x}_1}}=f\cdot {\displaystyle \frac{s_{Sp}}{\cos \left({\beta }_{I1}\right)}}\cdot I_C(x_1)+\left(1-f\right)\cdot {\displaystyle \frac{s_{Sp}}{\cos \left({\beta }_{I2}\right)}}\cdot J_C(x_1)-\\ -2\left(a_{Sp}+s_{Sp}\cdot \left(1-f\right)\right)\cdot I_S(x_1)+2s_{Sp}\cdot \left(1-f\right)\cdot J_S(x_1)\end{array} \\ \begin{array}{l}{\displaystyle \frac{d{J}_S(x_1)}{d{x}_1}}=-\left(1-f\right)\cdot {\displaystyle \frac{s_{Sp}}{\cos \left({\beta }_{I1}\right)}}\cdot I_C(x_1)-f\cdot {\displaystyle \frac{s_{Sp}}{\cos \left({\beta }_{I2}\right)}}\cdot J_C(x_1)-\\ -2s_{Sp}\cdot \left(1-f\right)\cdot I_S(x_1)+2\left(a_{Sp}+s_{Sp}\cdot \left(1-f\right)\right)\cdot J_S(x_1)\end{array} \end{gathered}$$

(39)

where $a_{Sp}$ and $s_{Sp}$ are light absorption and scattering coefficients in the spongy mesophyll.

System of Equation (39) was transformed to the system of Equation (40) to simplify analysis:

$$\begin{gathered} {\displaystyle \frac{d{I}_C(x_1)}{d{x}_1}}=L_{11}\cdot I_C(x_1) \\ {\displaystyle \frac{d{J}_C(x_1)}{d{x}_1}}=L_{22}\cdot J_C(x_1) \\ {\displaystyle \frac{d{I}_S(x_1)}{d{x}_1}}=L_{31}\cdot I_C(x_1)+L_{32}\cdot J_C(x_1)+L_{33}\cdot I_S(x_1)+L_{34}\cdot J_S(x_1) \\ {\displaystyle \frac{d{J}_S(x_1)}{d{x}_1}}=L_{41}\cdot I_C(x_1)+L_{42}\cdot J_C(x_1)+L_{43}\cdot I_S(x_1)+L_{44}\cdot J_S(x_1) \end{gathered}$$

(40)

where coefficients (L) correspond to coefficients of light flows in the system of Equation (39).

We used the method of undetermined coefficients describing light flows as the combination of elementary exponents ($I_C(x)=A{e}^{\lambda x_1}$, $J_C(x)=B{e}^{\lambda x_1}$, $I_S(x)=C{e}^{\lambda x_1}$, and $J_S(x)=D{e}^{\lambda x_1}$)

Equation (41) is the characteristic equation of this system:

$$\left(L_{11}-\lambda \right)\cdot \left(L_{22}-\lambda \right)\cdot \left[\left(L_{33}-\lambda \right)\cdot \left(L_{44}-\lambda \right)-L_{34}{L}_{43}\right]=0$$

(41)

The solution of this equation is Equation (42):

$$\begin{gathered} {\lambda }_1=L_{11} \\ {\lambda }_2=L_{22} \\ {\lambda }_3={\displaystyle \frac{L_{33}+L_{44}}{2}}+\sqrt{{\left({\displaystyle \frac{L_{33}+L_{44}}{2}}\right)}^2-\left(L_{33}\cdot L_{44}-L_{34}\cdot L_{43}\right)} \\ {\lambda }_4={\displaystyle \frac{L_{33}+L_{44}}{2}}-\sqrt{{\left({\displaystyle \frac{L_{33}+L_{44}}{2}}\right)}^2-\left(L_{33}\cdot L_{44}-L_{34}\cdot L_{43}\right)} \end{gathered}$$

(42)

Thus, the system of Equation (43) was used to describe the optical properties of the spongy mesophyll layer (for x ranging from h to h + l):

$$\begin{gathered} I_C(x)=A_1{e}^{{\lambda }_1\left(x-h\right)}+A_2{e}^{{\lambda }_2\left(x-h\right)}+A_3{e}^{{\lambda }_3\left(x-h\right)}+A_4{e}^{{\lambda }_4\left(x-h\right)} \\ {J}_C(x)=B_1{e}^{{\lambda }_1\left(x-h\right)}+B_2{e}^{{\lambda }_2\left(x-h\right)}+B_3{e}^{{\lambda }_3\left(x-h\right)}+B_4{e}^{{\lambda }_4\left(x-h\right)} \\ {I}_S(x)=C_1{e}^{{\lambda }_1\left(x-h\right)}+C_2{e}^{{\lambda }_2\left(x-h\right)}+C_3{e}^{{\lambda }_3\left(x-h\right)}+C_4{e}^{{\lambda }_4\left(x-h\right)} \\ {J}_S(x)=D_1{e}^{{\lambda }_1\left(x-h\right)}+D_2{e}^{{\lambda }_2\left(x-h\right)}+D_3{e}^{{\lambda }_3\left(x-h\right)}+D_4{e}^{{\lambda }_4\left(x-h\right)} \end{gathered}$$

(43)

where A₁, A₂, A₃, A₄, B₁, B₂, B₃, B₄, C₁, C₂, C₃, C₄, D₁, D₂, D₃, and D₄ are constants.

System of Equation (44) was used for the description of boundary conditions:

$$\begin{gathered} A_1+A_2+A_3+A_4=I_C(h) \\ {B}_1{e}^{{\lambda }_{1}l}+B_2{e}^{{\lambda }_{2}l}+B_3{e}^{{\lambda }_{3}l}+B_4{e}^{{\lambda }_{4}l}=J_C(h+l) \\ {C}_1+C_2+C_3+C_4=I_S(h) \\ {D}_1{e}^{{\lambda }_{1}l}+D_2{e}^{{\lambda }_{2}l}+D_3{e}^{{\lambda }_{3}l}+D_4{e}^{{\lambda }_{4}l}=J_S(h+l) \end{gathered}$$

(44)

I_C(x) and J_C(x) cannot be dependent on other light flows. It means that A₂ = A₃ = A₄ = 0 and B₁ = B₃ = B₄ = 0; in contrast, A₁ = I_C(h) and B₂ = J_C(h + l).

Based on elementary exponents and Equation (40) for I_S(x) and J_S(x), the system of Equation (45) was derived:

$$\begin{gathered} C\cdot \left(L_{33}-\lambda \right)+D\cdot L_{34}=-\left(A\cdot L_{31}+B\cdot L_{32}\right) \\ C\cdot L_{43}+D\cdot \left(L_{44}-\lambda \right)=-\left(A\cdot L_{41}+B\cdot L_{42}\right) \end{gathered}$$

(45)

Equations (46) and (47) are solutions of this system at λ₁ (A₁ = I_C(h) and B₁ = 0):

$$C_1=I_C(h)\cdot {\displaystyle \frac{-L_{31}\cdot \left(L_{44}-{\lambda }_1\right)+L_{41}\cdot L_{34}}{\left(L_{33}-{\lambda }_1\right)\cdot \left(L_{44}-{\lambda }_1\right)-L_{34}\cdot L_{43}}}$$

(46)

$$D_1=I_C(h)\cdot {\displaystyle \frac{-L_{41}\cdot \left(L_{33}-{\lambda }_1\right)+L_{31}\cdot L_{43}}{\left(L_{33}-{\lambda }_1\right)\cdot \left(L_{44}-{\lambda }_1\right)-L_{34}\cdot L_{43}}}$$

(47)

Equations (48) and (49) are solutions of this system at λ₂ (A₂ = 0 and $B_2=J_C(h+l)\cdot e^{-{\lambda }_{2}l}$):

$$C_2=J_C(h+l)\cdot e^{-{\lambda }_{2}l}\cdot {\displaystyle \frac{-L_{32}\cdot \left(L_{44}-{\lambda }_2\right)+L_{42}\cdot L_{34}}{\left(L_{33}-{\lambda }_2\right)\cdot \left(L_{44}-{\lambda }_2\right)-L_{34}\cdot L_{43}}}$$

(48)

$$D_2=J_C(h+l)\cdot e^{-{\lambda }_{2}l}\cdot {\displaystyle \frac{-L_{42}\cdot \left(L_{33}-{\lambda }_2\right)+L_{32}\cdot L_{43}}{\left(L_{33}-{\lambda }_2\right)\cdot \left(L_{44}-{\lambda }_2\right)-L_{34}\cdot L_{43}}}$$

(49)

Equations (50) and (51) describe relations between C₃ and D₃ for λ₃ and C₄ and D₄ for λ₄:

$$D_3=-C_3\cdot {\displaystyle \frac{L_{33}-{\lambda }_3}{L_{34}}}$$

(50)

$$D_4=-C_4\cdot {\displaystyle \frac{L_{33}-{\lambda }_4}{L_{34}}}$$

(51)

Based on the system of Equation (44), the system of Equation (52) was derived:

$$\begin{gathered} C_3+C_4=I_S(h)-C_1-C_2 \\ {D}_3{e}^{{\lambda }_{3}l}+D_4{e}^{{\lambda }_{4}l}=J_S(h+l)-D_1{e}^{{\lambda }_{1}l}-D_2{e}^{{\lambda }_{2}l} \end{gathered}$$

(52)

where C₁, C₂, D₁, and D₂ can be calculated with Equations (46)–(49). Combining Equations (50)–(52), we derived the system of Equation (53):

$$\begin{gathered} C_3+C_4=I_S(h)-C_1-C_2 \\ {C}_3{\displaystyle \frac{L_{33}-{\lambda }_3}{L_{34}}}{e}^{{\lambda }_{3}l}+C_4{\displaystyle \frac{L_{33}-{\lambda }_4}{L_{34}}}{e}^{{\lambda }_{4}l}=-J_S(h+l)+D_1{e}^{{\lambda }_{1}l}+D_2{e}^{{\lambda }_{2}l} \end{gathered}$$

(53)

Equations (54) and (55) are solutions of this system:

$$C_3={\displaystyle \frac{\left(I_S(h)-C_1-C_2\right)\cdot \left(L_{33}-{\lambda }_4\right)e^{{\lambda }_{4}l}-\left(-J_S(h+l)+D_1{e}^{{\lambda }_{1}l}+D_2{e}^{{\lambda }_{2}l}\right)\cdot L_{34}}{\left(L_{33}-{\lambda }_4\right)e^{{\lambda }_{4}l}-\left(L_{33}-{\lambda }_3\right)e^{{\lambda }_{3}l}}}$$

(54)

$$C_4={\displaystyle \frac{\left(-J_S(h+l)+D_1{e}^{{\lambda }_{1}l}+D_2{e}^{{\lambda }_{2}l}\right)\cdot L_{34}-\left(I_S(h)-C_1-C_2\right)\cdot \left(L_{33}-{\lambda }_3\right)e^{{\lambda }_{3}l}}{\left(L_{33}-{\lambda }_4\right)e^{{\lambda }_{4}l}-\left(L_{33}-{\lambda }_3\right)e^{{\lambda }_{3}l}}}$$

(55)

Combining Equations (50), (51), (54) and (55), we derived Equations (56) and (57) to calculate D₃ and D₄:

$$D_3=-{\displaystyle \frac{L_{33}-{\lambda }_3}{L_{34}}}{\displaystyle \frac{\left(I_S(h)-C_1-C_2\right)\cdot \left(L_{33}-{\lambda }_4\right)e^{{\lambda }_{4}l}-\left(-J_S(h+l)+D_1{e}^{{\lambda }_{1}l}+D_2{e}^{{\lambda }_{2}l}\right)\cdot L_{34}}{\left(L_{33}-{\lambda }_4\right)e^{{\lambda }_{4}l}-\left(L_{33}-{\lambda }_3\right)e^{{\lambda }_{3}l}}}$$

(56)

$$D_4=-{\displaystyle \frac{L_{33}-{\lambda }_4}{L_{34}}}{\displaystyle \frac{\left(-J_S(h+l)+D_1{e}^{{\lambda }_{1}l}+D_2{e}^{{\lambda }_{2}l}\right)\cdot L_{34}-\left(I_S(h)-C_1-C_2\right)\cdot \left(L_{33}-{\lambda }_3\right)e^{{\lambda }_{3}l}}{\left(L_{33}-{\lambda }_4\right)e^{{\lambda }_{4}l}-\left(L_{33}-{\lambda }_3\right)e^{{\lambda }_{3}l}}}$$

(57)

2.5. Description of Several Iterations of the Light Propagation Through Leaf

Equations from Section 2.2, Section 2.3 and Section 2.4 can be used for the first iteration of the calculation of light transmittance and reflectance. However, Equations (13), (14), (25) and (26) show that light can secondarily input into the leaf lamina; that is, the second iteration of calculation of the light propagation is possible (as well as the third iteration, fourth iterations, fifth iteration, etc.). This effect can be large at low light absorption coefficients (particularly for the near-infrared light, NIR [3]).

The second iteration of calculation of the light propagation through the leaf and light transmittance from lamina to air was based on the same equations from Section 2.2, Section 2.3 and Section 2.4, which were used for the first iteration, after substitution of the following parameters:$I_C(0)={I_C}^{RI}$, $J_C(h+l)={J_C}^{RI}$, $I_S(0)={I_S}^{RI}$, $J_S(h+l)={J_S}^{RI}$, ${\beta }_{I1}={\beta }_{I2}$, and ${\beta }_{I2}={\beta }_{I1}$.

Calculated after this procedure I_C^RI, J_C^RI, I_S^RI, and J_S^RI could be used for the third iteration of calculation, etc. Thus, Equations (58) and (59) could be used for the calculation of light outputs from adaxial and abaxial leaf surfaces (J_out and I_out, respectively):

$$J_{out}={\displaystyle \sum _{i=1}^N{J_{out}}^i}$$

(58)

$$I_{out}={\displaystyle \sum _{i=1}^N{I_{out}}^i}$$

(59)

where J_outⁱ and I_outⁱ are light outputs from adaxial and abaxial leaf surfaces, respectively, which are calculated on iteration i, and N is the quantity of the iterations, which is necessary to approximately describe the reflectance and transmittance spectra of leaves.

2.6. Description of Light Absorption Coefficients

In the current model, we considered that the light absorption coefficient was the function of the light wavelength through photosynthetic pigments and complexes formed by these pigments [49,56,57,58]. Maier et al. [49] assumed that the concentration of photosynthetic pigments in the spongy mesophyll layer was 20% of this concentration in the palisade mesophyll layer (N_Sp/P = 0.2):

$$a_{Sp}\left(\lambda \right)=N_{Sp/P}\cdot a_P\left(\lambda \right)$$

(60)

where $a_P\left(\lambda \right)$ and $a_{Sp}\left(\lambda \right)$ are spectra of the light absorption coefficient of layers of the palisade and spongy mesophyll, respectively. $a_P\left(\lambda \right)$ and $a_{Sp}\left(\lambda \right)$ were used as $a_P$ and $a_{Sp}$ in the analysis.

Equation (61) was used for the calculation of $a_P\left(\lambda \right)$:

$$a_P\left(\lambda \right)=C_{ChA}\cdot a_{ChA}\left(\lambda \right)+C_{ChB}\cdot a_{ChB}\left(\lambda \right)+C_{Car}\cdot a_{Car}\left(\lambda \right)$$

(61)

where C_ChA, C_ChB, and C_Car are concentrations (mg cm⁻³) of chlorophyll a, chlorophyll b, and carotenoids, respectively, in the palisade mesophyll, and $a_{ChA}\left(\lambda \right)$, $a_{ChB}\left(\lambda \right)$, and $a_{Car}\left(\lambda \right)$ are spectra of specific light absorption coefficients (cm² mg⁻¹) of chlorophyll a, chlorophyll b, and carotenoids, respectively. These spectra were constructed on the basis of [58] (Figure 3a).

Equation (62) was used for the calculation of the average concentration of each photosynthetic pigment (C^av) in leaves to compare with experimental concentrations (Figure 3b):

$$C^{av}=C\cdot \left({\displaystyle \frac{h}{l+h}}+N_{Sp/P}\cdot {\displaystyle \frac{l}{l+h}}\right)$$

(62)

where C is the pigment concentration in the palisade mesophyll layer.

2.7. Experimental Methods Used to Parameterization and Verification of the Model

Peas (Pisum sativum L., variety “Albumen”), which are dicot plants, were used to parameterize and verify the developed model of leaf reflectance. Plants were cultivated in pots with universal soil (nine plants per pot) under conditions of a vegetation room (16 h photoperiod and 24 °C); plants were irrigated three times a week. Luminescent lamps FSL YZ18RR (Foshan Electrical And Lighting Co., Ltd., Foshan, China) were used for the illumination. There were two variants of the illumination intensity (Figure S1) as follows: moderate (about 55 µmol m⁻²s⁻¹) and low (about 15 µmol m⁻²s⁻¹); these intensities were controlled using the Thorlabs PM100D optical power meter with an S120VC sensor (200–1100 nm) (Thorlabs Inc., Newton, MA, USA). Here, 9-day-old and 16-day-old pea plants were investigated in different variants of the experiment. Moreover, 16-day-old pea plants, which were cultivated under the 55 µmol m⁻²s⁻¹ light intensity, were used to parameterize the model. All experimental variants were used to verify the developed model.

Measuring leaf reflectance spectra, content of chlorophylls and carotenoids, and thickness of the palisade and spongy mesophyll layers were analyzed. Second mature leaves were investigated. All noted measurements were carried out on the same leaves.

Reflectance spectra were measured in attached pea leaves before other measurements. The handheld PolyPen RP 410 UVIS system (Photon Systems Instruments, Drásov, Czech Republic), which is the specialized system to measure spectral characteristics of leaf reflectance in plants, was used.

Cross-sections of pea leaves were manually prepared using a razor blade in accordance with the widely used botanical method. The leaf fragment was fixed by two plates of foam plastic; further, cross-sections were cut, placed into a water drop on the microscope slide, and fixed by the cover glass. Microscope MT 4200 with lens Plan Achromat 40X and ocular micrometer (Meiji Techno, Saitama, Japan) was used to estimate thicknesses of the palisade and spongy mesophyll layers.

Concentrations of chlorophyll a, chlorophyll b, and carotenoids were measured using the standard biochemical method [59,60,61] of spectrophotometry of acetone extract from pea leaves using an SF-2000 UV/Vis Spectrophotometer (OKB Spectr, St. Petersburg, Russia). Optical densities were measured at 440.5, 644, and 662 nm; concentrations of photosynthetic pigments were calculated on the basis of Holm–Wettstein equations.

3. Results

3.1. Parameterization of the Developed Leaf Reflectance Model

Parameterization of the developed model of leaf reflectance was the first task. Thicknesses of palisade and spongy mesophyll layers (Figure 2b) and concentrations of the chlorophyll a, chlorophyll b, and carotenoids (Figure 3b) were estimated using experimental methods in the leaves of pea plants (see Section 2.7). Measured reflectance spectra of these leaves were used to check the accuracy of the developed model. It should be noted that pea leaves were used as an example of leaves of dicot plants, which have palisade and spongy mesophyll layers [38]; the model cannot be used to describe reflectance of leaves of monocot plants with uniform mesophyll.

The following parameters were used for initial parameterization. Average experimental values of thicknesses (h = 35.5 μm and l = 58.6 μm) were used. C_ChA = 2.77 mg cm⁻³, C_ChB = 1.69 mg cm⁻³, and C_Car = 0.94 mg cm⁻³ were used because average concentrations of chlorophyll a, chlorophyll b, and carotenoids in leaves, which were calculated on basis of these C_ChA, C_ChB, and C_Car, Equation (62), and N_Sp/P = 0.2 (47), corresponded to experimental ones. $a_{ChA}\left(\lambda \right)$, $a_{ChB}\left(\lambda \right)$, and $a_{Car}\left(\lambda \right)$ were constructed on the basis of [58] (Figure 3a). Values of other parameters were β_O1 = 35º (in accordance with the angle of the leaf illumination performed using PolyPen RP 410), I₀ = 1000 μmol m⁻²s⁻¹ and J₀ = 0 μmol m⁻²s⁻¹ (assumed), F_S = 0 (assumed), n_I = 1.415 [48], f = 0.5 (assumed), =5 cm⁻¹, and =1000 cm⁻¹ [47].

Figure 4 shows model-based reflectance spectra, which were calculated at different quantities of the iterations (N) in accordance with Equation (58). It was shown that reflectance of visible light was accurately described at N = 2 (or more) because increasing N did not influence the reflectance spectra at N ≥ 2. Moreover, N = 1 could be used for the calculation of the leaf reflectance in blue and red spectral regions; a small error was observed in the green spectral region only. In contrast, an accurate description of the NIR reflectance was observed at N ≥ 5 because increasing N influenced model-based leaf reflectance spectra at N ≤ 5. As a result, we used N = 6 for further analysis to completely exclude technical error, which could disrupt results at low N.

Figure 5 shows experimental and model-based spectra of leaf reflectance. The model-based spectrum was calculated at initial values of parameters (see above) and N = 6. This result showed that the model-based spectrum was similar to the experimental one in a qualitative manner; however, quantity differences were also observed (especially in blue, green, and NIR spectral regions).

An additional parameterization of the developed model was the next task. First, we analyzed the influence of $s_{Sp}$ on the accuracy of description of experimental reflectance spectra (Figure 6). It was shown that decreasing $s_{Sp}$ from 1000 cm⁻¹ to 600 cm⁻¹ increased this accuracy (Figure 6a–c); in contrast, decreasing s_Sp from 600 cm⁻¹ to 400 cm⁻¹ decreased similarity between model-based and experimental spectra (Figure 6c). As a result, we used =600 cm⁻¹ in further analysis.

Second, we used small corrections of C_ChA and C_ChB to increase the determination coefficient between experimental and model-based reflectance spectra. It was shown that using C_ChA = 3.19 mg cm⁻³ and C_ChB = 2.09 mg cm⁻³ improved the accuracy of the description of the experimental leaf reflectance spectrum (Figure 7a). Average leaf concentrations of chlorophyll a and b, which were calculated on the basis of these values of C_ChA and C_ChB, were similar to experimental ones (differences approximately equaled to standard errors of experimental values of concentrations of chlorophylls). Thus, we used corrected C_ChA and C_ChB in further analysis of the model.

Third, we analyzed the influence of the fraction of the rough surface (F_S) on reflectance spectra to decrease the error of the model-based leaf reflectance spectra in the blue spectral region. It was shown (Figure 8) that increasing F_S to 0.1–0.225 improved description of the experimental leaf reflectance spectrum (especially in the blue spectral region). As a result, we used F_S = 0.15 in the further analysis of the model.

Thus, additional parameterization showed that relatively minor changes in parameters provided increasing the determination coefficient between experimental and model-based leaf reflectance spectra from about 0.93 to about 0.99. This improvement was especially large in blue, green, and NIR spectral regions.

3.2. Verification and Analysis of the Developed Leaf Reflectance Model

The next task of the investigation was the verification of the developed model on the basis of experimental leaf reflectance spectra. We used plants, which were cultivated under moderate (about 55 µmol m⁻²s⁻¹) and low (about 15 µmol m⁻²s⁻¹) light intensities, to provide varying leaf characteristics because light intensity is an important factor influencing leaf anatomy and the content of photosynthetic pigments [62,63,64]. Two ages of pea plants (9 and 16 days) were used to additional check results.

It was shown (Figure 9) that cultivation under low light intensity significantly decreased thicknesses of palisade and spongy mesophyll layers and increased the content of carotenoids in 16-day-old pea plants. The concentrations of chlorophylls were not dependent on the light intensity during cultivation. Figure 10 shows experimental and model-based leaf reflectance spectra. It was experimentally shown that cultivation of pea plants under low light intensity increased reflectance in the green spectral range and decreased reflectance in the NIR spectral range (Figure 10a). Analysis of the developed model showed similar results as follows: low thicknesses of palisade and spongy mesophyll layers (equaling to experimental thicknesses in plants cultivated under low light intensity) increased reflectance in the green spectral range and decreased reflectance in the NIR spectral range (Figure 10b). In contrast, changes in the carotenoid concentration weakly influenced the model-based leaf reflectance spectrum (Figure S2). The last result showed that the difference in the carotenoid concentration was not a probable reason for changes in the leaf reflectance spectrum in pea plants cultivated under low light intensity.

A similar experimental investigation of 9-day-old pea plants showed (Figure S3) that cultivation under the low light intensity did not influence concentrations of chlorophylls and carotenoids; in contrast, thicknesses of palisade and spongy mesophyll layers were decreased in this case. Decreased thicknesses of palisade and spongy mesophyll layers were accompanied by increasing reflectance in the green spectral range and decreasing reflectance in the NIR spectral range in both experimental and model-based reflectance spectra (Figure S4). This result additionally supported the key role of changes in thicknesses in changes of reflectance because the carotenoid concentration was not influenced by cultivation under the low light intensity in the last experiment.

Thus, it was shown that the developed model could simulate the dependence of the leaf reflectance spectrum on thicknesses of palisade and spongy mesophyll layers. The similarity of changes was observed for two experimental groups (9-day-old and 16-day-old pea plants). This result verified the developed leaf reflectance model and showed the possibility of using the model to analyze factors influencing reflectance.

It should be finally noted that changes in the model-based leaf reflectance spectrum were shown at simultaneous decreasing thicknesses of palisade and spongy mesophyll layers (in accordance with the experimental decreasing). Next, we used the developed model to analyze the influence of changes in thickness of the palisade mesophyll layer only and in thickness in the spongy mesophyll layer only. It was shown (Figure 11a) that decreasing the thickness of the palisade mesophyll layer increased reflectance in the green spectral region but did not influence reflectance in the NIR spectral region. In contrast, decreasing the thickness of the sponge mesophyll layer both increased reflectance in the green spectral region and decreased reflectance in the NIR spectral region (Figure 11b).

Thus, the last results of the analysis of the developed model showed that decreasing thicknesses of both palisade and spongy mesophyll layers stimulated reflectance in the green spectral region. In contrast, decreasing the reflectance in the NIR spectral region was caused by decreasing the thickness of the spongy mesophyll layer only.

4. Discussion

Multispectral and hyperspectral imaging is an important tool of remote sensing of plants that can be used to estimate their growth and development [12,13,14,15], photosynthesis [16,17,18,19,20], nitrogen content [11,12], water content [8,10], concentrations of pigments [5,6,7], and other characteristics, including revealing changes in these parameters under the action of stressors. Potentially, using reflectance in narrow bands and calculating reflectance indices should increase efficiency of the plant remote sensing [1,3]; however, relationships between reflectance in specific narrow spectral bands (and, therefore, corresponding reflectance indices) and specific plant characteristics can be unstable (e.g., strong variability was shown for relationships of photochemical reflectance index to photosynthetic parameters [22,23]). Possible reasons for the instability include changes in photosynthetic pigment composition, variations in leaf anatomy, and different angles between the direction of the incident light and the leaf surface [3].

The development of mathematical models describing the optical properties of leaves provides theoretical tools to reveal the influence of noted reasons on reflectance spectra and to minimize this influence. The main result of the current study is the development, parameterization, and verification of the analytical model to describe leaf reflectance spectra in dicot plants. This model describes light flows through the leaf surface on the basis of Snell’s and Fresnel’s laws [51,53,65]. Fractions of smooth and rough surfaces, which provide transmittance and reflectance collimated and scattered light [65,66], are also considered in the developed model. Light flows in the palisade mesophyll, which has a low light scattering coefficient [49], are mainly described on the basis of the Beer–Bouguer–Lambert law. Light flows in the sponge mesophyll, which has a high light scattering [49], are described on basis of the Kubelka–Munk theory (in accordance with [40,42]). Parameterization and verification of the developed model (with using pea leaves) show that this model accurately describes leaf reflectance spectra (see Figure 8c, Figure 10 and Figure S4), that is, it can be used to analyze data of plant remote sensing on the basis of multispectral and hyperspectral imaging.

The developed model does not require to assume that leaves are described as series of plates; this assumption is the basis of numerous plant optical models, including different variants of PROSPECT and FLUOSPECT [38,66,67,68]. Our model does not require a complex description of light trajectories, which is typically for ray tracing models [38,55,69]. Moreover, combining different approaches (particularly the Beer–Bouguer–Lambert law for the palisade mesophyll and the Kubelka–Munk theory for the spongy mesophyll), the developed model provides a relatively simple description of the optical properties of leaf; the proposed equations have the analytical solution. It should be additionally noted that the developed model can be easily adapted to describe reflectance in monocot plants with one type of leaf mesophyll (h = 0 µm and N_Sp/P = 1 should be used).

The additional result of the current study is the theoretical revealing of factors that influence the leaf reflectance in different spectral ranges. Particularly, reflectance in the NIR spectral range can be stimulated by increasing the light scattering in the spongy mesophyll (Figure 6) and thickness of this mesophyll (Figure 10 and Figure S4). These effects are in good accordance [29] with the positive influence of the ratio of mesophyll cell surface to intercellular air spaces on reflectance in the NIR spectral region because increasing this ratio should stimulate light scattering in the spongy mesophyll layer; increasing the spongy mesophyll thickness can also increase the NIR reflectance. Influence of the light scattering and thickness of the spongy mesophyll layer can cause variability of vegetation reflectance indices, which use reflectance in the NIR spectral range, including the normalized difference vegetation index [70], optimized soil-adjusted vegetation index [71], green normalized difference vegetation index [72], triangular vegetation index [73], and others.

Reflectance in the green spectral range can be mainly stimulated by increasing the light scattering in the spongy mesophyll (Figure 6) and decreasing thicknesses of both the palisade and spongy mesophyll layers (Figure 10 and Figure S4). These effects can increase the variability of indices that use the green spectral bands, including the photochemical reflectance index [16,17,18,19,20] and its modifications [74,75,76]. It is probable that the increased variability of indices can explain relatively weak relationships between the photochemical reflectance index and photosynthetic parameters [22,23].

Finally, reflectance in the blue spectral region is relatively stable; however, increasing the fraction of the rough surface stimulates reflectance in this region (Figure 8). Considering this effect, it can be expected that reflectance indices, which are used reflectance in the blue spectral region (e.g., normalized difference pigment index [77], normalized phaeophytinization index [78], Carter index 1 [79], and others), are dependent on ratio between areas with smooth and rough leaf surfaces. It should be additionally noted that the leaf surface reflectance is mainly considered as the Fresnel’s reflectance [65]; however, results of parameterization of our model (on basis of pea leaves) show that a fraction of the rough leaf surface can be about 15% and more (Figure 8c,d).

Thus, the developed model of leaf reflectance in dicot plants shows some potential reasons for the variability of reflectance indices and, therefore, the variability of their relationships to plant characteristics. In future investigations, the model can be used both for the complex analysis of the leaf reflectance spectra and for the revealing of new reflectance indices, which will be stably related to plant characteristics.

It should be finally noted that the current model has some assumptions and restrictions that should be considered at its use. (i) The model described two layers of mesophyll (palisade and spongy), which are typical for dicot plants [38]. It means that equations of the model can be potentially used to analyze reflectance of leaves of different dicot plants; however, the model parameterization is mainly based on experimental results (contents of photosynthetic pigments, thicknesses of layers of palisade and spongy mesophyll, and leaf reflectance spectra) shown on pea plants. Thus, quantity analysis of the leaf spectra in other dicot plants can require additional parameterization of the developed model.

(ii) The model does not consider participation of the epidermal layer in the light reflectance and transmittance because this layer is narrow and transparent (see, e.g., Figure 2a); it can be expected that it weakly influences light propagation. However, epidermal cells can also influence leaf optical properties because they focus light on mesophyll cells and, probably, induce additional light scattering [80,81]. The description of the effect is one of the potential ways of model development in the future.

(iii) It is known that the low light intensity is accompanied by stomatal closure [82] and chlorophyll content changes [83,84] modifying photosynthetic processes. Other environmental factors, for example, temperature or humidity, can also influence stomata opening, content of photosynthetic pigments, and photosynthetic parameters [82,85,86]. The current model is focused on the leaf optical properties and does not describe these physiological processes (Figure 1). However, considering the relationship of plant reflectance spectra to concentrations of photosynthetic pigments [3,4,7], photosynthetic activity [17,19,22,23], and water content [8,9], it is probable that description of changes in concentrations of pigments, photosynthesis, and transpiration under action of environmental factors can be integrated to the leaf reflectance model in the future; this “extended” model can provide an effective tool for analysis of results of the plant remote sensing in the future.

5. Conclusions

Improving methods of plant remote sensing on the basis of multispectral and hyperspectral imaging requires the development of mathematical models of reflectance in plants. In the current study, we developed the analytical model of the leaf reflectance in dicot plants (on example of pea leaves), which was based on the Snell’s and Fresnel’s laws to describe light reflectance and transmittance on borders of lamina, on the Beer–Bouguer–Lambert law to describe the light transmittance in the palisade mesophyll, and on the Kubelka–Munk theory to describe light transmittance and scattering in the spongy mesophyll. The model was parameterized and verified using the experimental results and the literature data. Analysis of the model theoretically showed some factors influencing the leaf reflectance spectra (e.g., the coefficient of light scattering in the spongy mesophyll, thicknesses of both mesophyll layers, fraction of the rough surface in the leaf, and others).

In future investigations, the developed model can be used both for the complex analysis of leaf reflectance spectra to improve interpretation of results of remote sensing and for revealing new reflectance indices, which will be stably related to plant characteristics.