Uniﬁed Four-Stream Radiative Transfer Theory in the Optical-Thermal Domain with Consideration of Fluorescence for Multi-Layer Vegetation Canopies

: Vegetation radiative transfer models (RTMs) are important tools to understand biosphere-atmosphere interactions. The four-stream theory has been successfully applied to solve the radiative transfer problems in homogeneous canopies for both incident solar radiation, thermal and ﬂuorescence emission since 1984. In this note, we describe the development of a uniﬁed radiative transfer theory for optical scattering, thermal and ﬂuorescence emission in multi-layer vegetation canopy, and provide a detailed mathematical derivation for the ﬂuxes inside and leaving the canopy. This theory can be used to develop vegetation models for remote sensing applications and plant physiological processes, such as photosynthesis and transpiration. It can also be used to solve the radiative transfer problems in soil-water, soil-water-atmosphere, or soil-vegetation-atmosphere ensembles, besides the soil-vegetation system presented in the note.


Introduction
Vegetation radiative transfer models (RTMs) describe the interaction between light and vegetation canopies, and they simulate radiative fluxes within and above vegetation canopies. The radiation inside vegetation canopies interacting with leaves is a prime driver and control of many important plant physiological processes, such as photosynthesis and transpiration. These processes determine the exchange of energy and matter between vegetation and the atmosphere. Hence, RTMs are an indispensable component for modeling and understanding plant physiology, and the hydrology and micro-meteorology of soil-vegetation-atmosphere systems. The leaving radiation over vegetation canopies can be observed by using remote sensing techniques, and the intensities and spectral characteristics of the observed signals are affected by leaf biochemical and biophysical properties, canopy structure, and sun-observer geometry. Therefore, RTMs connecting the vegetation characteristics with remote sensing observations can be applied to support the quantitative use of remote sensing data for vegetation monitoring.
A number of canopy RTMs have been developed for studies of energy budget of vegetated ground surface and for the interpretation of remote sensing data. Classical canopy RTMs simulate the directional scattering and absorption of radiation by canopy elements like leaves or needles and compute the canopy radiation budget and outgoing radiation as an outcome of these scattering and absorption events. These models require spectral signatures of leaves, which are either measured or calculated by a leaf sub-model. The PROSPECT model by Jacquemoud and Baret [1] is the most common one, but other options are available. Geometric optics (GO) models are another group of RTMs which put more emphasis on describing vegetation canopy radiative transfer by means of simplified opaque sub-canopies of various geometric shapes (e.g., cubic, cylindrical, conical, and ellipsoidal) and the spatial distributions of the sub-canopies [2]. They simulate radiation leaving canopies or land surfaces with probabilities of seeing components in a scene, such as sunlit and shaded crowns, and sunlit and shaded soil. Radiances of different scene components may be estimated empirically, with separate classical RTMs or given as input. The combination of GO and classical RTMs results in GORT (Geometric optics and radiative transfer) hybrid models [3]. GO models are widely used for remote sensing applications. However, for supporting the simulation of plant physiological processes, the classical RTMs or the GORT models are preferred, as the GO models usually pay much less attention to the scattering and absorption of radiation by leaves within vegetation canopies, where most physiological processes occur.
Vegetation RTMs differ from each other mainly by representation of the canopy and solutions of radiative transfer problems. Canopies generally exhibit large horizontal and vertical heterogeneity of both leaf biophysical properties and canopy structural properties as adaptations to environmental conditions. The large or even infinite variability of these properties requires simplifications in modeling the radiative transfer. A vegetation canopy can be represented simply by a turbid medium [4][5][6], a multilayer medium [7][8][9], or explicitly by a realistic 3-D natural scene [10,11]. Canopy radiative transfer problems can be solved using computer simulation techniques, such as Monte Carlo ray tracing [11,12] or analytically using techniques, like the discrete-ordinate method [13]. These have resulted a number of vegetation RTMs, such as the Suits model [5], the SAIL (Scattering by Arbitrarily Inclined Leaves) model [6], the DART (Discrete Anisotropic Radiative Transfer) model [10,14], and the FLIGHT (Forest Light Interaction) model [11].
The four-stream theory, as a special case of the discrete-ordinate method, is the most popular scheme being used in radiative transfer modeling [6,15,16]. It accounts the interactions of four types of fluxes with a vegetation canopy, namely upward and downward diffuse isotropic fluxes, a direct solar flux, and a flux associated with the radiance in the direction of observation. The four-stream theory allows modeling the directional variations of diffuse fluxes, and can be applied to leaves, canopies, atmosphere and water bodies. For example, both the Suits and SAIL models employ the four-stream theory and a turbid medium representation of vegetation canopies. The SAIL model has been successfully applied in remote sensing in the optical domain, in which the only radiation source is the solar radiation. Later, it has been extended to thermal domain [17] and adapted for fluorescence radiative transfer [18][19][20]. As a result, a unified four-stream theory has been established to solve radiative transfer problems of optical radiation, thermal and fluorescence emission for homogeneous canopies. Models based on the unified theory have obvious benefits in interpreting multi-source remote sensing data, since they share a common architectural description.
Furthermore, the Soil Canopy Observation of Photosynthesis and Energy fluxes (SCOPE) model incorporates this unified theory and provides predictions of radiation within vegetation canopy, which is linked with several processes sub-models [19]. For example, the net photosynthetically active radiation (PAR) on individual leaves, as an output of radiative transfer modeling, is used as input of the FvCB (Farquhar-von Caemmerer-Berry) photosynthesis model to simulate leaf photosynthetic production [21]. However, the above mentioned applications (Suits, SAIL, and SCOPE) assume canopies as homogeneous turbid media. To investigate the effects of vertical heterogeneity of leaf properties on reflectance, photosynthesis and fluorescence, Yang et al. [9] adapted the four-stream theory to multi-layer vegetation canopies for the radiative transfer of optical radiation and fluorescence emission and and developed the mSCOPE model (i.e., multi-layer SCOPE). In Yang et al. [9], new separated solutions to the radiative transfer problems in the optical domain, and for emitted fluorescence, are presented for multi-layer canopies, while, for the radiative transfer in the thermal domain, the solution in SCOPE is kept, which involves an iteration approach. After further investigating, we find that one solution can be applied to the radiative transfer problems in the optical-thermal domain from 0.4 µm to 50 µm, as well as for fluorescence emission from 640 nm to 850 nm. This led to a unified radiative transfer theory.
In the present note, we focus on the explicit derivation of the four-stream radiative transfer theory in multi-layer vegetation canopies when thermal and fluorescence emission of radiation by leaves and the soil is incorporated. Besides the solution for multi-layer canopies in Section 3, we also describe the applications of the four-stream theory in homogeneous canopies in Section 2. Both sections start with subsections that summarize of the general idea and itemized derivations in Sections 2.1 and 3.1. In Section 3.4, we present some illustrative simulations to show the capability of the unified theory. Comparison with field measurements is not included in this note but can be found in the work of mSCOPE [9]. We expect this note could help the readers to better understand the use of the four-stream theory in vegetation radiative transfer modeling and the radiative transfer theory in SAIL, SCOPE, and mSCOPE.

Generalized Form of the Four-Stream Theory without Emission
The four-stream theory considers four types of fluxes involved in the radiative transfer processes, namely a direct solar flux, two hemispherical (semi-isotropic) diffuse fluxes (up-and downward) and a flux in the direction of viewing. Their corresponding (equivalent) irradiance is E s for direct solar flux, E − for downward diffuse flux, E + for upward diffuse flux, and E o for flux in observer direction.
Canopies are represented by a horizontally homogeneous and continuous canopy layer with gaps. Hence, the vertical variation of fluxes is considered and computed while the fluxes in the horizontal plane are constant except for the hot spot effect. Following the Kubelka-Munk theory, the four-stream theory expresses the radiative transfer in the vertical direction as a set of linear differential equations [6]. Each differential equation describes changes in irradiance (E) of the four fluxes with respect to a change in the height (x) due to absorption and scattering. They are given in matrix-vector notation as: where x is the vertical relative height to the canopy bottom surface (i.e., top of canopy: x = 0; bottom of canopy: x = −1), and L is canopy leaf area index (LAI). The coefficients in the matrix are the extinction and scattering coefficients. The extinction coefficients (k and K) are determined by gap probabilities which depend on sun-leaf and observer-leaf geometry, respectively. The geometries can be defined with the zenith angles (θ) and azimuth angles (ϕ) of the sun, leaf, and observer. The scattering coefficients (s, a, σ, s , w, v, v ) depend not only on these geometries but also the optical characteristics of foliar elements, namely leaf reflectance ρ and transmittance τ. Note that the main symbols used are listed at Appendix as 'List of Symbols'. Furthermore, vegetation canopies are usually bounded with a soil layer beneath. The upward outgoing fluxes of the soil surface include reflection of incident fluxes on the soil, which is the flux emanating from the canopy at the bottom: where R xx comprises the soil surface reflectance factors. The two-letter subscripts indicate the incident and outgoing fluxes types: d refers to the diffuse fluxes, s to the direct solar flux, and o to the flux in The first two steps refer to the properties of optically thin vegetation layers with LAI of Ldx and for an isolated thick vegetation layer with LAI of L (i.e., the entire canopy layer), respectively. The third step is to obtain the optical properties of a combined system of soil and the thick vegetation layer. The optical properties of either a vegetation layer or a soil-vegetation system do not require solving the fluxes explicitly at all levels and can be obtained based on the general solutions of the differential equations. For example, the development of the Suits and SAIL models covers the first three step and provides canopy reflectance simulation of a soil-vegetation system. However, advanced applications of the four-stream theory, such as photosynthesis and energy balance simulations in the SCOPE model, require the fluxes in the vertical direction E(x) or even the radiation on individual leaves. This requires the particular solution of differential equations mentioned in step 4, which can be achieved by employing the boundary conditions: incident direct solar irradiance E sun and incident diffuse sky irradiance E sky .

Scattering and Extinction Coefficients
The first step for solving the radiative transfer problems is to compute scattering and extinction coefficients based on the inputs of the angles of the sun (θ s , ϕ s ), leaves (θ l , ϕ l ), and observer (θ o , ϕ o ), and leaf reflectance and transmittance. To compute these coefficients of a thin layer (Ldx), it is essential to simulate the interaction of the four fluxes with the leaves of various orientations.
Leaf and layer area projections to the directions of the leaf normal, sun and observer are necessary inputs to translate radiometric quantities of individual leaves into those the of a thin layer. With these projections, Verhoef [6] gave the scattering efficiency factors Q sc (E 1 , E 2 ), where E 1 refers to the type of incident flux, and E 2 to the type of scattered or attenuated fluxes. Here, we summarize the expressions for the coefficients, and for the detailed derivation the reader is referred to Reference [6]. Tables 1 and 2 present the symbols and expressions, respectively, of the scattering and extinction coefficients for possible combinations of incident (E 1 ) and scattered or attenuated (E 2 ) fluxes of single leaves. Table 1. Symbols of scattering and extinction coefficients for the possible combinations of incident (E 1 ) and scattered or attenuated (E 2 ) fluxes of single leaves.
to the viewing direction. f o = cos θ l + tan θ o sin θ l cos(ϕ l − ψ). f 1 and f 2 are the leaf area projections to layers, where f 1 refers to the greater of the two. f 1 = (1 + cos θ l )/2 and f 2 = (1 − cos θ l )/2 . ρ and τ are reflectance and transmittance of leaves, respectively.

Scattering Matrix of an Vegetation Layer
The early applications of the four-stream theory are mainly for reflectance simulation without generating any output in the form of diffuse fluxes inside the canopy [6,15]. For reflectance simulation, it is essential to find the scattering matrix (reflectances and transmittances) of the vegetation layer that links fluxes at top of canopy with those at bottom of canopy: where b and t refer to bottom and top of the layer (i.e., b referring to x = −1, and t referring to x = 0), respectively. τ xx and ρ xx are the reflectance and transmittance factors of the layer. The two-letter subscripts indicate the incident and outgoing flux types.
To obtain the scattering matrix, one has to employ the general solutions of the differential equations. In a vertically homogeneous canopy, both the leaf optical properties and canopy structural characteristics (e.g., leaf angles) are vertically identical. Note that, in the horizontal plane of thin layers (Ldx), there might be leaves with different orientations quantified by a leaf inclination distribution function (LIDF), but the LIDF remains unchanged in the vertical direction. As a result of these assumptions, the coefficients in Equation (1) are independent of vertical relative height (x). This leads to an analytical solution for the linear differential equations. For example, the general solution for where C 1 is a constant. Therefore, we can obtain τ ss = E s (−1)/E s (0) = exp(−kL). Symmetrically, τ oo = exp(−KL). The hemispherical reflectance and transmittance properties of the isolated canopy layer are derived by linking the diffuse fluxes at the lower and upper boundaries of the layer. A complete derivation of the scattering matrices is given in Reference [15].

Reflectance of a Soil-Vegetation System
Once the scattering matrix is computed, the form of four-stream theory in the soil-vegetation system can be obtained by combining Equation (3) and Equation (2) and is given as: where and The flux interaction with a soil-vegetation system is summarized in Equation (5)  To solve the system in Equation (5) for the effective reflectance at top of the layer (or the surface reflectance of the system) R * t , we introduce Equation (5b) into Equation (5a) and obtain that where R * t describes the relationship between the downward and upward fluxes at a level while taking into account the multiple reflections with the surface under the layer of interest: i.e., E u (t) = R * t E d (t).
Here, we introduced a new quantity T * d , which may be called the effective downward transmittance, since it describes the relationship between the downward fluxes at two successive levels while taking into account the multiple reflections with the surface under the layer of interest: . Therefore, we obtain two important quantities for a surface-layer system that link the fluxes at top and bottom of the layer.
Top: These two quantities show the background surface leads to greater reflectance and transmittance of the system than the isolated layer itself (i.e., R * t > r t and T * d > t d ). For 'black background' (i.e., R b = 0), the effective reflectance and downward transmittance become the reflectance and transmittance of the vegetation layer: R * t = r t and T * d = t d .

Extension of Four-Stream Theory to the Thermal Domain
Verhoef et al. [17] extended the four-stream theory to the thermal domain by adding an extra source from emitted thermal fluxes by leaves and soil.
where ε v is the single-leaf emissivity. The hemispherical flux H v (x) is the result of thermal emission by blackbody leaves at a temperature T v and is given by where B is Planck's spectral radiance function for the indicated leaf temperature and the associated wavelength. Assuming all the leaves in the canopy have the same temperature and emissivity, the extra terms of Equation (8) are independent of vertical height (x) and thus the differential equations (Equation (8)) can be solved in the same way as the classic four-stream differential equations (Equation (1)). However, in natural canopies temperature may vary vertically, and a different solution to the radiative transfer problem should be applied, which is discussed in section 3.
Note that the sum of leaf emissivity (ε v ), transmittance (τ) and reflectance (ρ) at the same wavelength is unity according to Kirchhoff's law of thermal radiation. Hence, leaf emissivity can be expressed by using leaf reflectance and transmittance. Similar to Equation (3), the connection between fluxes at top of canopy and bottom of layer including the thermal emission, can be given as follows: Solving the differential equations (Equation (8)), we can obtain the hemispherical emissivity (γ d ) and effective emissivity (γ o ) in the viewing direction of the isolated canopy layer, which can be found in Reference [22]. Furthermore, comparing the emissivity with the reflectance and transmittance of vegetation layer, we find that According to Equation (11), γ d and γ o can also be interpreted as the hemispherical and directional absorptance, respectively, since the sum of reflectance, transmittance and absorptance is unity. Therefore, this obeys Kirchhoff's law of thermal radiation for entire vegetation layer, which is a proof of the derivation.
For the soil-vegetation system, Equation (5) is revised as: where ε s is the soil emissivity. The hemispherical flux H s is the result of thermal emission at soil temperature T s . Introducing Equation (12b) into Equation (12a), we obtain the surface reflectance R * t , as well as emitted thermal radiation, at top of the canopy U t (t).
The system shown as Equation (12) (both (12a) and (12b)) is the general form of Equation (5) and can be applied in the entire optical-thermal domain. In the solar-reflective, optical regime, the emitted fluxes from leaves (H v ) and soil (H s ) are often neglected. Hence, Equation (12) is simplified to Equation (5), or, vice versa, Equation (5) is a special case of Equation (12). This suggests that the four-stream theory can be unified for the optical and thermal domain for homogeneous canopies, namely by using Equation (12).

Four-Stream Theory in Multi-Layer Canopies
The radiative transfer problem of vertically heterogeneous canopies is more complex. Although the general form in Equation (1) still holds for multi-layer canopies, the vertical heterogeneity of leaf properties and leaf orientations as illustrated in Figure 2b may lead to vertical variation of the scattering and extinction coefficients. As a result, the differential equations cannot be solved directly as done for vertically homogeneous canopies. Therefore, the analytical expression of the scattering matrix (Equation (3)) is difficult to derive. Furthermore, both fluorescence and thermal emission can vary vertically due to the varying illumination. Consequently, the approach proposed in Reference [17] does not hold in this case, and both the system shown as Equation (12) and the corresponding solution should be revised. The radiative transfer in multi-layer canopies for top-of-canopy (TOC) reflectance simulation is given in Reference [15] by using the adding method. Yang et al. [9] extended the adding method for calculating the radiative flux profiles in the canopy, which allows simulation of fluorescence and photosynthesis. In what follows, we provide a unified theory for the radiative transfer in both the optical and thermal domain with the contribution of fluorescence.

Four-Stream Theory in Multi-Layer Canopies with Fluorescence and Thermal Emission
Any vertically heterogeneous vegetation layer can be regarded as a thick layer consisting of N thin homogeneous layers. The value of N should be sufficiently large to ensure each thin layer is vertically homogeneous and the LAI of one thin layer is small enough (i.e., δL < 0.1). Therefore, a practical value of N could be 10 times canopy LAI. We use the term 'elementary layers' for these thin layers.
In a heterogeneous N-layer system that is bounded by a surface at the bottom, we may distinguish the N layers by a number from 1 to N, and the fluxes at the bottom and the top of the system by the numbers 1 and N + 1. The levels at the interfaces between layers are numbered from 2 to N (Figure 2c). Using this numbering for layers and their interfaces, the following set of equations describes radiative transfer in the whole system: . . .
The radiative transfer of fluxes can be summarized in a similar way as in Figure 1, but now we consider that each layer may generate fluxes (E u e and E d e ) and the fluxes at the top of a surface may contain a reflected downward flux, as well as an upward emitted flux, called U. Figure 3 is a revised version of Figure 1 with contribution of thermal or fluorescence emission.
The following quantities are required to solve the system in Equation (14)  Because the computation of fluorescence emission normally requires knowing the excitation fluxes inside the canopy, the radiative transfer of fluorescence can only be solved when the radiative transfer of incident radiation is solved (i.e., E u and E d in the third item). Likewise, the temperatures are sometimes modeled as a function of the optical fluxes inside the fluxes (i.e., E u and E d in the third item), for example, with an energy balance approach as in the SCOPE model [19]. Nevertheless, the leaf temperatures are sometimes given as direct input. In this case, the radiative transfer of canopy emitted thermal radiation can be combined with radiative transfer of incident radiation.
The radiative transfer for either the incident, emitted fluorescence, and thermal radiation can all be represented by the system expressed as Equation (14). If only incident radiation is considered, then the emission of each layer and soil is set to zero. If the fluorescence and thermal emission from soil and foliage are zero (i.e., E d e , E u e and U(1) are zero), then it follows that the upward emitted fluxes at every level U(j) all equal to zero. Therefore, the above presented theory is a unified radiative transfer theory for optical, thermal, and fluorescence in multi-layer canopies. In what follows, we first describe the general solution of such a system assuming thermal and fluorescence emission is known and then elaborate the computation of thermal and fluorescence emission afterwards.

Scattering Matrices of an Elementary Layer
For a thin layer, the scattering matrices can be directly computed from the extinction and scattering coefficients.
where δL is the LAI of a thin layer, which is canopy L/N. Note that all the scattering and extinction coefficients can vary vertically, as well as the reflectances and transmittances of the layer in the scattering matrix.

Solutions for a Surface-Layer System with Emission
The differences of fluxes at any two successive levels are affected by the thin layer between the two levels, and the thick layer under the lower level. This can be regarded as a surface-layer combination with surface reflectance R * t (j) and upward emitted fluxes U(j) and the layer scattering matrices (t d (j), r b (j), r t (j), t u (j)), and the fluxes interactions between two successive levels can be expressed as: The relationship between the downward fluxes at successive levels, while taking into account both the multiple reflectance with and the emission from the thick layer, can be obtained by introducing Equation (16b) into Equation (16a): where we find For the emission in the viewing direction, we employ the coefficient K that projects a small leaf area on the plane perpendicular to the viewing direction. Therefore, the conversion of the diffuse emitted radiation E b e and E f e of individual leaves into those for the layer can be expressed as where ε l and T l are the leaf emissivity and temperature, respectively. E + et and E − et are upward and downward diffuse irradiance, and E et o is equivalent thermal irradiance in the viewing direction. E et s is direct solar thermal irradiance emitted from the leaves, which is null for leaves and set to zero in the calculation. It is given in the above equations for the sake of the competence of the four-stream theory.
The upward emitted thermal fluxes by soil surface can be computed in the same way as where ε s and T s are the soil emissivity and temperature, respectively.

Fluorescence Emission
Fluorescence emission of a thin layer is excited by the direct solar flux (E s ), upward (E + ), and downward (E − ) diffuse light at each spectral band (λ e ) from 400 to 700 nm. One can simulate the upward and downward fluorescence emission for a given irradiance of a thin layer and a fluorescence emission efficiency. The source terms due to fluorescence emission by leaves are given as where the excitation-emission coefficients (with subscript f ) are determined by sun-observer geometry, canopy structure, leaf optical properties, and fluorescence emission efficiency of photosystems, as given in Reference [9]. The fluxes (E s , E + and E − ) at each elementary layer are computed as described in Section 3.2. The excitation-emission coefficients with subscript f can be found in Table 3, which, by analogy, are related to the corresponding scattering coefficients in Table 2. The excitation-emission matrices M b and M f are analogous to leaf reflectance and transmittance, respectively. The difference is that reflectance and transmittance describe the relationship between incident radiation and outgoing scattered radiation of a leaf, while the matrices describe the relationship between incident radiation and outgoing emitted radiation of a leaf. In a scattered event, the wavelength of an incident radiation and outgoing radiation remains the same, but, in a emission event, the emitted radiation is of different wavelength than the incident. Table 3. Emission coefficients for the possible combinations of incident (E 1 ) and emitted fluorescence (E 2 , with notation "ef") flux of single leaves.
and M f are the excitation-emission matrices of leaves in the back and forward direction (relative to the incident light), respectively.

Illustrative Simulations
To demonstrate the capability of the present theory, we carry out an illustrative simulation. The setting of the scenarios is briefly described as follows. A vegetation canopy with two distinctive layers, of which the upper layer consists of leaves with high chlorophyll content (C ab = 80 µg cm −2 ) and the lower with C ab = 20 µg cm −2 . The emissivity of the leaves in the canopy is set to 0.95. Moreover, the temperature of the leaves in the upper layer is set to 25 • C and the lower 15 • C. Besides leaf chlorophyll content and temperature, other properties of the two layers are identical. Specifically, both layers have LAI of 1.5 (i.e., total canopy is 3) and have spherical leaf angle distribution. The sun-observer geometry is defined as θ s = 30 • and θ v = 0 • .
The unified four-stream theory resolves the radiative transfer in vegetation canopies, and the profiles for fluorescence and thermal emission, upward, and downward diffuse fluxes, and vegetation absorption, as well as outgoing directional radiation, are simulated as shown in Figure 4. Different leaf chlorophyll content and leaf temperature of the two layers result in different leaf optical properties ( Figure 4a) and thermal emission (Figure 4b), respectively. Furthermore, the fluorescence emission clearly shows vertical variation (the spectral region from 640 to 850 nm in Figure 4b). This is mainly affected by vertical variation of the radiation inside the vegetation canopy as shown in Figure 4c,d. The radiation in the photosynthetically active radiation (PAR) region for the downward diffuse flux decreases from top to bottom due to the strong absorption in this region (Figure 4c). The absorption of PAR decreases gradually from the top to the bottom of the canopy (Figure 4e). An abrupt change of absorbed PAR (APAR) occurs in the middle of the canopy (i.e., x = 1/2) because of the change of high chlorophyll content of the upper layer to low chlorophyll content of the lower layer. Lastly, the radiance observed in the viewing direction is simulated, as well (Figure 4f), which can be used to compute canopy reflectance. Note the simulation results cover the spectral region from 0.4 µm to 50 µm and are in 1 nm spectral resolution, but other spectral resolutions are also possible for this unified theory.

Hot Spot Effects for Remote Sensing Applications
This note presents a theory for the simulation of all the four streams, including fluxes in the viewing direction, which supports remote sensing applications. The fluxes inside the canopy can be used to predict net radiation and photosynthesis. Hence, this theory can be applied as a component in vegetation energy balance, photosynthesis, or physiological models.
However, in this report, the hot spot effect on fluxes in the viewing direction (E o ) is not considered. The hot spot effect is caused by the preferential reflectance (or emission) of a surface in the opposite direction to the incident radiation. This leads to an enhanced brightness in that direction relative to other directions. In a vegetated surface consisting of optically isotropic leaves and soil, the hotpot is caused by the coincidence of exposure to incident radiation and the direct transmittance of scattered radiation in opposite direction. In other words, if the incident radiation reaches a leaf by direct transmittance, then the scattered radiation can directly escape the canopy following the reversed path. The probability of direct transmittance may also be enhanced in directions close to this path, i.e., angles close to the incident angle. The range of angles in which the direct transmittance is enhanced determines the size of the hotspot (in angular units), which increases with the size of the objects (leaves) relative to the vegetation height. The radiative transfer theory presented in this report uses a turbid medium representation of the vegetation. In this representation, leaves are infinitely small, and the preferential reflectance only occurs the exact incident direction, such that the hotspot is infinitely small.
Accounting for the variation of E o in the horizontal plane caused by the hot spot effects in vegetation with leaves of a determined size requires further geometrical modeling and accounting for co-variances of projections in incident and viewing direction. One can compute E o separately by combining the four-stream theory with Kuusk's approach for the hot spot effect [16,23].

Vertical Distribution of Leaf Area
We have considered the vertical variation of leaf optical properties and leaf orientations, which are expressed by the scattering and extinction coefficients of layers. In contrast, the distribution of leaf density in the vertical direction was not considered. Figure 2 depicts a uniform vertical distribution of leaves in the canopy, but a non-uniform arrangement of leaves in the vertical direction is very common in reality.
The radiative transfer theory presented here uses units of leaf area index (Ldx), rather than height above the soil surface, to define the depth of the canopy. The solution of the differential equations (Equation (1)) is therefore indifferent to the ways that the leaf area is distributed in the vertical direction. Thus, the use of a turbid medium representation of the canopy leads to the counterintuitive conclusion that the radiative transfer is indifferent to the vertical arrangement of leaf density.
However, one can imagine cases in which this is not the case, and the vertical distribution of leaf area matters for the radiative transfer. First, if optical properties are assigned to the air space between the leaves, then the optical properties of a layer are no longer determined by the properties of leaves alone. In that case, the vertical distribution of leaf density will affect the properties of the layers because the contribution of air to the layer properties will scale inversely with the leaf area density. Since the theory allows for layer properties to vary vertically, this can be accounted for easily by assigning appropriate layer properties. Second, if the objective is to model the hotspot or the position of sunlit and shaded locations in the canopy, then the turbid medium representation does not suffice. In that case, the theory can be blended with an approach to model specific effects of geometrical shapes. Apart from radiative transfer, the vertical distribution of leaf area also affects the physiology of vegetation via its effect on turbulent heat and matter exchanges and hydraulic pressure in the vegetation.

Applications to Other Systems
The theory is an extension of the classical four-stream theory to accommodate vertical heterogeneity of leaf properties and orientations in vegetation canopies (i.e., soil-vegetation systems). The modeling concept, however, has potential in radiative transfer modeling in media other than vegetation, such as water and atmosphere, which have a clear multi-layer structure [16,24]. The main differences are in the derivation of scattering, extinction, and excitation-emission coefficients of the sub-layers. Once these coefficients are derived for each layer, the solution of the system can be solved by using the solutions presented in this report (Equations (17) and (18)).
Additionally, the unified four-stream theory can be used to solve the radiative transfer problems in soil-water, soil-vegetation, soil-water-atmosphere, or soil-vegetation-atmosphere ensembles [25]. These ensembles can be considered as a multi-layer system. For example, compared with the soil-vegetation system presented in this report, the application of the unified theory to a soil-water system in coastal areas can be achieved by replacing the vegetation layer with a water layer. The radiative transfer problems in a soil-vegetation-atmosphere can be solved by adding an atmosphere layer on top of a soil-vegetation system. Author Contributions: Conceptualization, software, writing-original draft preparation, writing-review and editing, P.Y.; conceptualization and writing-review and editing, W.V.; conceptualization, funding acquisition and writing-review and editing, C.v.d.T. All authors have read and agreed to the published version of the manuscript.

Funding:
This research was funded by the Netherlands Organization for Scientific Research, grant ALWGO.2017.018.

Conflicts of Interest:
The authors declare no conflict of interest.The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Symbol
Description Remarks E irradiance, radiant flux density spectral  Table 1