A Four-Scale Bidirectional Reflectance Model Based on Canopy Architecture 

Abstract

1  Introduction

In Li and Strahler's BRD models [22], a forest stand is assumed to consist of randomly-distributed objects containing leaves as turbid media. These two-scale models mark a major advancement in simulating radiation regimes in forest stands as compared with the one-scale turbid-media models, but also dramatically increase the complexity of the models. However, two-scale models may be regarded as much simplified mathematical descriptions of the physical reality. First, trees are generally not randomly distributed in space but are usually clustered at large scales due to variations in the soil and topographic conditions, creating patchiness of forest stands. They are also not randomly positioned when close to each other because of the natural repulsion effect in competition for resources (Franklin et al. [11]). Second, leaves are not randomly distributed within tree crowns. In conifer stands, for example, needles are grouped into shoots, branches and whirls, and all these sub-canopy structures are important in determining the radiation regime, especially the directional reflectance. In the present paper, these two additional scales of plant canopy architecture are considered. Figure  1 shows the concept of the different radiation transfer models based on canopy architecture at several scales. The objective of this paper is to present the first four-scale BRD model for the purpose of investigating the effect of the different architectural scales on the directional reflection behaviour of plant canopies. 
 

2  Model description

2.1  Tree distribution

According to the conditional probability theory, the probability of having i trees given j groups in a quadrat, P(i|j), times the probability of having j groups in the quadrat, P(j), gives the probability of having i trees in the quadrat given j groups. The probability of having i trees in the quadrat, P(i), is the summation of all the conditional probabilities for j:

P(i) = ¥ å j P(i|j)P(j).

(2)

When both P(j) and P(i|j) are determined by the Poisson process, we obtain

PN(i;m1;m2) = e-m1 mi2 i! ¥ å j = 1 [ m1 e-m2 ]j j! ·ji  for i = 0, 1, 2, ... ,

(3)

where 

m1 =  m2 n- m ,

(4)

is the mean number of groups per quadrat and

m2 =  n-m m .

(5)

is the cluster mean size. m = m₁ m₂ is the mean and n = m₁ m₂[1+m₂] is the variance of the distribution of the number of trees per quadrat and can be measured in a real forest stand. In creating this distribution, the following implicit assumptions are used [13]: Within the modelling domain each quadrat is equally likely to receive a cluster and the placement of a cluster is independent of the placement of any other clusters. 

1) The parameter m₁ corresponds to a priori value of the density of a cluster and m₂ corresponds to a priori value of the mean size of clusters (optional). 
2) In conformity with a Poisson model the variance about the mean size of cluster is equal to the mean size of cluster. 
3) The assignment of a cluster size to one location is independent of the assignment of any other cluster size or cluster location (optional). 
4) The points in a cluster are propagated by a `progenitor' who is located at the site of each cluster (pseudo-contagious assumption).

If the real distribution is not known, m₂ will be fixed at a low value ( m₂ = 1 implies a tree distribution close to the random case). Comparison with field data will later determine the values for those parameters. 

Figure 3 shows a random tree distribution (Poisson) in comparison with Neyman distributions with the mean grouping m₂ of 1, 5, and 20 trees per group. In the calculation, the average density was 75 trees per quadrat (1500 trees/ha with 20 quadrats of 500 m² each). Groupings of 1 and 5 produce distributions centred at 75 and increase the standard deviation from the random case. m₂ = 1 means that the grouping is mostly random (group of one) but it contains probabilities of having groups of two trees or more and therefore has a broader distribution than the random case. If the group size is too large compared to the overall mean (m), the Neyman distribution becomes variable with peaks at multiples of that size. As a result, the curve for m₂ = 20 in Fig. 3 shows maxima at 20, 40, and 60 trees per quadrat. At i = 0, the number of quadrats (having no trees) increases with group size, suggesting that for larger group sizes, the probability of having empty quadrats is greater. In Fig. 2, measurements from a jack pine forest are compared with the Poisson and Neyman distributions. The domain studied was one hectare with 1793 trees divided into 100 quadrats. A Neyman distribution with grouping of 3 simulates the measurements more closely than the Poisson theory. Processes other than Neyman could be used in this kind of analysis. For its simplicity and for the first analysis of tree clumping in forest canopies, only the Neyman, i.e., the double-Poisson process is used in this paper. 

2.2  Tree crown projection

A cone casts a shadow on the ground in addition to the cylinder beneath it only if the solar zenith angle is larger than one-half of the apex angle. This shadow area is expressed as [19]:

Sgc =  é ê ë 1 tan(gs) - p 2 +gs ù ú û r2 ,

(6)

where g_s (see Fig. 5 ) is defined as

gs = sin-1 æ ç è tan(a) tan(qs) ö ÷ ø ,

(7)

where q_s is the solar zenith angle. The shadow cast by the cylinder on the ground has an area given by

2tan(qs)Hb r .

(8)

To obtain the total area of a tree projected on the ground surface, the base of the tree, a disc, should also be included:

pr2 .

(9)

In this model, the cylinders can either rest on the ground surface or have their bases elevated to the same height (H_a). Since the lower part of a forest is generally a trunk space without much foliage, it is more realistic to elevate the bases, i.e. to have the crown ``on a stick''. In our model, the height of tree crown base affects the contribution of the ground surface to the hotspot. The summation of (6), (8) and (9) gives the total area within the outline of a tree crown projected on the ground surface :

Sg =  ì ï ï í ï ï î pr2 + 2tan(qs)Hbr for       qs£a,  é ê ë 1 tan(gs) + p 2 +gs ù ú û r2 + 2tan(qs)Hb r  for       qs > a.

(10)

This projected area can be used to find the ground area not seen by a viewer, if the viewer is far enough so that the parallax won't change the view. For this reason, we only need to replace q_s in (10) with q_v. This quantity is denoted by V_G, which is the area on the ground covered by opaque crowns in the viewer's direction:

Vg =  ì ï ï í ï ï î pr2 + 2tan(qv)Hbr  for       qv£a, é ê ë 1 tan(gv) + p 2 +gv ù ú û r2 + 2tan(qv)Hb r  for       qv > a.

(11)

After calculating the shadow area that a single tree casts on the ground, the sunlit crown proportion seen by the viewer is computed from the total surface area of the tree visible to the viewer projected to a plane perpendicular to the view line. For the cylinder we have

tab = 2rsin(qv)Hb .

(12)

The total tree surface area projected on a plane perpendicular to the view line is 

tac =  ì ï í ï î pr2 for       qv = 0,  pr2cos(qv) for       qv < a,  pr2cos(qv)+tact for       qv > a.

(13)

t_act of (13) is described in Appendix A. 

For cylinders, the self-shadowing and illumination geometry are very simple: half the surface is illuminated and half in shadow. The illumination part seen by the viewer is then

tib = 2rsin(qv)Hb· é ê ë 1- f p ù ú û ,

(14)

where f is the azimuth angle difference between the sun and the viewer. The sunlit part of a cone seen by a viewer, denoted t_ic, cannot be easily expressed in a simple equation, but is described in Appendix A. This quantity depends on the apex angle of the trees and the positions of the sun and the viewer. If q_s is smaller than a, no self-shadow occurs on the cone, but as q_s gets larger than a, a portion of the cone is shaded. The proportion of the tree that is illuminated and seen by the viewer (P_ti) can be found by adding the areas illuminated and dividing by the total surface area of the tree presented to the viewer. However, the contributions from both the cone and the cylinder part of a tree must be included. Thus

Pti =  tic + tib tac+tab .

(15)

2.3  Canopy gap fraction

The gap probability, P_gap(q), is well known for a continuous medium [28]. For a discontinuous canopy, we use an equation similar to Li and Strahler [21] but modified to consider the foliage clumping effect [4]:

Pgap(qs) = e-G(qs)Lo ·WE/gE ,

(22)

where 

Lo = m· _ s

(23)

and 

m = L/[V·D/B]

(24)

is the foliage density, and 

_ s (qs) =  V Sg(qs) cos(qs) ,

(25)

is the mean path length of the solar beam through a tree. V is the volume contained within the tree crown outline, D/B is the number of trees per m², W_Eis a clumping index of the shoots within tree crowns [4], L is the leaf area index (LAI) defined as half the total needle area per unit ground surface area [4]. 

(22) gives the probability of the solar beam passing through a single tree. Using q_v instead of q_s gives the probability of seeing the background through a single tree. For a canopy with a random foliage angle distribution, G(q) = 0.5. For canopies with branch architecture [3], G(q) = a-bq, where a and b are positive constants. Adding a gap probability to tree crowns gives less shadow on the ground and more ground area seen by the viewer. 

3  Mutual shadowing effect

Near the principal plane on the backscatter side, the illumination and observation of a sunlit area can occur in the same gap in the canopy. In such a case, the processes of illumination and viewing are correlated. The angle range over which the correlation occurs depends on the tree size and the average tree spacing. A function f(f), depending on the difference in the azimuth angle between the sun and the viewer, is used to consider this correlation effect. Using the mean distance between tree centres [] calculated by: 

E (d) =  1 2    ___ ÖD/B ,

(34)

we can determine the azimuthal angle range within which overlapping of one tree with the other occurs in both the illumination and view directions: 

tan(qE) =  2r E(d)-2r .

(35)

A linear equation is used between f = 0 and f = q_E: 

f(f) = 1 - f qE ,

(36)

where f(f) is a kernel with non zero values only when 0 < f < q_E. The overall mutual shadowing effects, Q_c and Q_b, for the conical and cylindrical parts, respectively, are then given by

Qc = Qvc ·Qsc + f(f)·[Qsvc-Qvc ·Qsc]

(37)

and

Qb = Qvb ·Qsb +f(f) ·[Qsvb-Qvb ·Qsb],

(38)

where Q_svc is the minimum of Q_sc and Q_vc, and Q_svb is the minimum of Q_sb and Q_sb. With considerations of these mutual shadowing effects, the probability of seeing illuminated tree crowns becomes

Pti =  Qc ·tic + Qb ·tib Qvc ·tac + Qvb ·tab .

(39)

Much of the complexity of the mutual shadow calculations arises from the fact that the conical and the cylindrical parts have different geometry and have to be separately modelled. Since conifer trees are better represented using the combination of these two parts than the single cone shape, we believe that the formulation presented above is an improvement over previous models based on cone geometry when used for conifer forests. 

4  Canopy gap size and hotspot

The theory of gap size distribution in plant canopies was first derived by Miller and Norman [25] and further investigated by Chen and Black [5] and Chen and Cihlar [6]. The accumulated gap size distribution can be calculated from:

Fat(l) =  é ê ë 1+Lt l Wt ù ú û e- Lt[1+[(l)/(Wt)]] ,

(40)

for gaps ³ l, where W_t is the characteristic width of a tree crown projected on the ground; and L_t is the clumping-adjusted projected tree crown area index. This formulation is made according to the finding of Chen and Cihlar [7] that a large part of the gap size distribution curve measured in boreal forests is determined by tree crowns as discrete objects. W_t and L_t are computed as follows: 

Wt =    __ ÖSg ,

(41)

and

Lt = Wt Sg D/B ,

(42)

where S_g is the area of a tree crown projected on the ground in the sun's direction and W_t is a tree clumping index, determined by the Neyman distribution as

Wt = log(Pig(Neyman))/log(Pig(Poisson))

(43)

It is 0.95 obtained from simulation with a Neyman grouping of 4 for a forest with 4000 stems per hectare divided into 25 quadrats.

Similar to (40), a distribution for gaps inside crowns can be computed where the shoots are taken as the foliage elements:

Fas(l) =  é ê ë 1+Ls l Ws ù ú û e- Ls[1+[(l)/(Ws)]] .

(44)

where

Ls =  G(qs)WE L gE cos(qs)

(45)

is the clumping-adjusted foliage area index and the characteristic width of the foliage elements is 

Ws =    __ ÖAs .

(46)

W_E is a clumping index for the effect of confining shoots within tree crowns; g_E is the needle-to-shoot area ratio quantifying the effect of needle grouping in shoots on the radiation interception; and A_s is the average shoot projected area. According to Chen [2], for the black spruce stand investigated in this paper, L = 4.5, G(q) = 0.5, W_E= 0.70, g_E = 1.41 and W_s = 35 mm. Only very small gaps are computed inside trees, but they contribute significantly to the total canopy gap fraction because of their large numbers. In the calculation of the hotspot, these small gaps are of critical importance in determining the shape of the hotspot.

[small modifications: equation (45) is use with LAI=1; W_Eis modified in the code since TRAC measurements include clumping between crowns. The clumping is reduced by half, e.g. 0.70 becomes 0.85]

4.1  Hotspot between crowns

4.2  Hotspot within tree crowns

The formulation presented above has implications on the inversion of the model for leaf area index because the micro-scale structures within the tree crowns affect the BRD and cannot be ignored.

In computing the hotspot within the tree crowns, the self-shadowing and the vertical tree crown structure must be considered. Discontinuities in mathematical expressions are found at nadir or the boundary between sunlit and shaded sides of a crown. Therefore, we treat the two sides separately. On the sunlit side, the probability of seeing the illuminated foliage P_T at the hotspot is simply 1 - P_ig, i.e. the viewer sees either the sunlit ground or the sunlit foliage. Considerable complication arises when we treat the tree crown surface as a complex surface with micro-scale structures. In this case, a view line can penetrate into the depth of the crown and reach to the shaded foliage behind the sunlit part even though the imaginary tree crown surface is theoretically sunlit. The problem may be understood in the extreme cases: 1) a tree crown is very dense - the imaginary surface can reasonably represent the tree crown, and 2) a tree crown is very sparse - the imaginary surface ceases to have meaning. The reality is in between these two extremes, i.e. the imaginary surface exists but is unsmooth or complex. In this case, the effect of mutual shadowing among foliage elements within a tree crown cannot be ignored.

Figure 10 shows the physics involved in the determination of the hotspot within a tree crown, where a tree is considered as a vertical structure with spherical foliage elements (shoots) dispersed within it. Using L(x) as the accumulated LAI from the imaginary surface to a given location x within the crown, the accumulated sunlit leaf area from 0 to L(x) is 1-e^{-C_s L(x)}, where C_s = G(q_s) W_E/g_E sin(q_s + a) to take into account the cone inclination and the vertical structure. At position x, the increment of sunlit leaf area with increasing L(x) is 

d dL(x) [1-e-Cs L(x)] = Cs e-Cs L(x)

(54)

Similarly, the increment of leaf area seen after dL(x) is

d dL(x) [1-e-Cv L(x)] = Cv e-Cv L(x).

(55)

If we denote L_H as the leaf area index accumulated horizontally from the sunlit side to the shaded side, the total probability of viewing sunlit leaf area from 0 to L_H is then

Q1(LH)  = G(x) ó õ LH 0 Cv Cse-CvL(x) e-CsL(x) dL(x)  = G(x)[1- e-LH(Cs +Cv)] é ê ë CsCv Cs +Cv ù ú û

(56)

where G(x) is the first-order scattering (geometric shadow) phase function of the foliage elements. It is defined as

G(x) =  é ê ë 1- Cp x p ù ú û

(57)

where C_p is a coefficient determined by the optical properties of foliage elements. If the elements are solid spheres with a Lambertian surface, C_p is unity, resulting in the phase function being 1.0 at x = 0, 0.5 at x = p/2 and 0 at x = p. The phase function in this case gives the proportion of the sunlit sphere surface seen by the viewer. Although a conifer shoot can be approximated by a sphere to describe the projected area [2], it can not be treated as a solid because of the gaps within it. The value of C_p for porous elements is smaller than unity. Because of lack of data, we assume that C_p = 0.75 in this paper. This constant is found to affect the sharpness of the hotspot. Q₁ must be calculated for all trees, it gives:

Q1 tot =  k å j = 1 Q1(j ·LH)·Ptj(Vg)

(58)

where P_tj(V_g) is the probability of a view line going through j tree crowns calculated by (20). 
[Q_1tot is no longer computed with (58). It uses the same reasonning as Q_2tot] On the shaded side, the probability of observing sunlit leaves also exists, especially in canopies with sparse tree crowns. The increment of sunlit leaf area along the sun's direction is calculated using (54), but the increment of leaf area viewed from the opposite side is calculated differently as follows:

d dL(x) [1-e-Cv(LH-L(x))] = -Cv e-Cv(LH-L(x)).

(59)

At a given depth x, the probability of observing sunlit foliage is

G(x)Cs e-Cs L(x)(-Cv e-Cv(LH- L(x))).

(60)

After integration from L_H to 0 across the tree crown, the overall probability becomes

Q2(LH) =  é ê ë 1- Cp x p ù ú û [e-CsLH- e-CvLH] é ê ë CsCv Cv - Cs ù ú û .

(61)

This phenomenon must be considered on each tree along the path of view. On the first tree, it is simply Q₂(L_H). For i trees along the path it is:

Q2 tot = Q2(LH)  k å i = 1 ·Ptt(i) ·K(i)

(62)

where

Ptt (i) = Pat(i) Pat(i-1) Pgapi-1(qv).

(63)

in which

Pat(i) =  k å j = i Ptj(Vg).

(64)

P_at(i) represents the probability of having i trees within one path. The probability of seeing the ith tree behind i-1 trees is found by the probability of having i-1 tree overlapping, P_at(i-1), times the probability of passing through i-1 trees, P_gap^i-1(q_v). As the view line penetrates through the forest, it reaches the lower portion of the canopy where shaded foliage is more likely to be observed. For the ith tree on the path, this height attenuation is considered by

K(i) = e-i GE(qv) ·cos(qv)/cos(qs).

(65)

With the proportion of sunlit and shaded tree crown as P_ti and (1-P_ti) respectively, we can then calculate the probability of seeing the illuminated foliage as follows:

PTf = Pti ·Q1 tot + [1-Pti]·Q2 tot.

(66)

This equation is only valid far from the hotspot. A hotspot function, similar to (51), can also be defined. F_s(x) follows the same principles as F_t(x) but it involves foliage elements (shoots) instead of tree crowns. 

Fs(x) =  ó õ ¥ lmins é ê ë 1- x tan-1(l/Hs) ù ú û Ns(l) dl ó õ ¥ lmins Ns(l) dl ,

(67)

where H_s is the effective distance between two layers of leaves inside a tree crown. In our model, it is inversely proportional to the leaf area density in the tree crown or L_o, i.e.

Hs =  r Lo ,

(68)

where r is the tree crown radius. This definition produces sharper hotspots in canopies with sparser tree crowns. After considering the hotspot on tree crowns, the overall probability of observing sunlit foliage is 

PT = PTf + [(1-Pig)-PTf]Fs(x) .

(69)

This gives P_Tf outside of the hotspot where F_s(x) = 0 and (1-P_ig) at the centre of the hotspot where F_s(x) = 1. Figure  9 shows examples of the hotspot kernel calculated using (51) and (67) along the principal solar plane at a solar zenith angle of 45^°. As the first approximation, the shape of the kernel is the same in all directions, i.e., the hotspot has a circular shape.

4.3  Canopy reflectance

5  Model results

5.1  Effect of Neyman grouping

The measured gap fraction data at large zenith angles are positively biased because of the effect of multiple scattering on the measurements using an optical instrument [2]. The modelled curves all show a sharp increase near the vertical direction because of the simple geometry to represent the tree crown and overlapping of the crowns without considering the repulsion effect. This creates unsmoothness of the model results at the nadir as shown later. In reality, the tree crown geometry is more variable and the reflectance distribution are usually smoother around nadir.

To understand the effect of the Neyman grouping, Fig. 12 shows different components of the model for a range of m₂ values from 1 to 50 at q_v = 30^°, q_s = 55^° and f = 0. The grouping of trees produces a larger gap fraction because it increases the probability of having quadrats with few trees. With the increasing grouping size, the proportion of canopy illuminated (P_ti) decreases, and so does the probability of seeing a tree illuminated (P_T). This is mainly because of the increase in the gap fraction. In accordance with P_vg, the probability of seeing the illuminated ground surface (P_G) increases. Red and near infra-red reflectances decrease with the Neyman grouping because ground reflectivity is low compared to the crown reflectivity.

5.2  BRD and hotspot

Comparisons of the model results were made with measurements of a boreal black spruce forest [9]. Figures 14a and 14b show the calculated and measured reflectance in the red band for two different q_s values, 40^° and 55^°, in the principal plane. At q_s = 40^°, the model is overestimating the reflectance on the forwardscattering side. This is mainly due to the constant reflectivities used for the shaded foliage and ground. The model uses constant values for F_dt and F_dg in (70) and (71) for the fraction of diffuse radiation in the total incident solar irradiance. As the fraction changes during the course of the day, these constant values better represent the average daily conditions. Near solar noon when the sun is high, the diffuse fraction is smaller than the daily average and the reflectivities for the shaded components are also correspondingly smaller than the daily average, resulting in the discrepancies between the modelled and observed values at q_s = 40^°. The comparisons suggest that the model will benefit from accurate separation of the diffuse and direct solar radiation. The model shows sharp spikes at the hotspot in comparison with the measurements. The gentle variation in the measurements may be a result of the low angular resolution (15^°) which effectively produced window-averaged results. The effect of the averaging is to dampen and broaden the peak at the hotspot. For q_s = 40^° and other small q_s values not shown here, the width of the hotspot is well modelled. The model does not perform as well at the larger solar zenith angles like Fig. 14(b) for q_s = 55^°. This model deficiency may have resulted from the simplified geometry used to represent the tree crowns. The model geometry is such that, all foliage is confined within the cone or cylinder, while in reality branches extend much further than the mean radius and intercept more radiation than the model prediction. One way to handle the problem is to increase the mean radius of the tree crowns, but in that case the branch structure needs to be more vigorously described to allow more gaps within the tree crown. This suggests that geometric-optical models are still approximations to reality and accurate simulation of BRD requires accurate descriptions of the canopy architecture at all scales. Figure 14c shows a comparison between the model and measurements for the reflectance in the near-infrared band at q_s = 40^°. The model is able to simulate the measurements closely. In Fig. 14d where q_s = 60^°, the model performs well except for the largest q_v, especially on the backscattering side, indicating the effect of non-uniform distribution of multiple scattering on BRD which is not considered in the present study.

Figure 15 shows hemispheric distributions of reflectance in the red band at four solar zenith angles. The hotspot varies in size, being generally larger when the sun is higher in the sky. The model does not show much of the usual bowl-shape because of the high LAI and the fixed diffuse fractions used in the model. Some of the bowl-shape can be seen in Fig. 15c and 15d. Figure 16 shows the corresponding distributions in the NIR band. The distributions are similar to those in Fig.  but the bowl-shape is more pronounced for q_s = 60^° and 75^°. The increase in reflectance at large q_v values, where the ground is hidden under the tree crowns is due in part to the larger reflectivity of the shaded foliage as compared to the shaded ground A singularity is seen at nadir for the large solar zenith angles, i.e. Fig. 16(d) because of the large shaded ground components viewed vertically.

6  Discussion

Much mathematical complexity arises, though still tractable, when we treat the tree crown surface as a complex surface within which mutual shadows and the hotspot also occur. Through our numerical simulation, we believe that the sub-canopy structures have profound effects on the directional reflection behaviour of the canopy and deserve such attention. Figure 17 shows the sensitivity of the model to the foliage density within tree crowns. With the fixed stand density and tree crown shape and dimensions (Table 1), the foliage density increases proportionally with LAI. For these simulations, we kept the multiscattering factors constant, and realised that some inaccuracies occurred in the model results because larger LAI should induce more multiscattering in the canopy. Both the red and NIR reflectances increase with decreasing LAI in the forward scattering direction due to the increased probabilities of observing sunlit foliage from the shaded side through sparse tree crowns. The hotspot is smaller for a sparser canopy (lower LAI) because of the increased probability of observing the ground, which has a lower reflectivity than the foliage. These model results support the findings of Soffer [32] that the measured reflectance on the shaded side of jack pine crowns is considerably larger than the predictions of a Li and Strahler model. Black spruce crowns are generally denser and the shaded side appears darker than the jack pine. Our model is able to simulate such differences.

The assumption of random tree distribution is removed at the expense of lengthy mathematical descriptions using the combined Neyman and binomial distributions. Apparently, such effort has not been very rewarding with regard to its impact on the final BRD results for the stand investigated. Figure 18a shows that the Neyman grouping increases the reflectance in the red band only on the forward scattering side near the vertical direction. The increased openness of the forest due to the Neyman grouping is the main cause of the increase in the modelled reflectance. The Neyman grouping decreases the width of the hotspot because a larger m₂ gives a smaller W_t which affects the hotspot kernel for the ground. At large view angles, a larger Neyman grouping gives a smaller reflectance because the contribution of the ground to the reflectance becomes very small and the tree crown properties dominate at these angles. The change of the hotspot width is more important in the NIR than the red band because the Neyman grouping affects the proportion of the ground illuminated (P_G). Our present simulation is limited to small quadrats with small Neyman groupings. More research is needed to determine the grouping effect at larger scales (large quadrats with large groupings).

Sub-canopy architecture is not only important in investigating the BRD of vegetated surfaces but also critical in the inversion of the model to obtain biophysical parameters such as leaf area index. Clumping of needles within shoots, for example, affects the radiation interception in the plant canopy and hence vegetation indices derived from the reflectances in red and near-infrared bands. Grouping of shoots in branches also has similar effects. Accurate estimation of LAI from the model requires appropriate descriptions of architecture at these sub-canopy levels.

In this paper, the effect of multiple scattering on BRD is not investigated. Our attention is first given to canopy architecture because the first order scattering is usually much larger than the sum of multiple scattering and the angular distribution of the multiply scattered radiances is usually considered to be isotropic [31]. However, at NIR wavelengths, the multiple scattering effect is larger than that at visible wavelengths and the isotropic assumption may lead to some inaccuracies in the BRD results. Such subtle differences may have significant effects on the simulated vegetation indices using reflectance factors for these two bands. We hold the believe that in modelling the multiple scattering, the canopy architecture is foremost important because it dictates the direction of first order scattering and the probability of observing the reflecting surfaces at different steps of the scattering sequence. The approaches of Goel et al. [14] and Strahler and Jupp [33] in the description of plant canopy architecture have such merits, but forest canopies are much more complex than what have been described. The model presented here provides a framework in which the effect of canopy architecture at various scales on BRD can be systematically investigated.
 

In our model, a tree crown consists of conical and cylindrical parts. An effective mathematical scheme is devised to estimate the mutual shadowing effect between tree crowns and among foliage elements within tree crowns. The canopy gap size distributions between and within tree crowns are used to describe the hotspot size and shape. The description, for the first time, closely relates the bidirectional reflection behaviour to canopy attributes.

8  Acknowledgements

List of symbols

A  Cone geometry

When the solar zenith angle q_s is less than the apex angle a, self-shadowing on the cone occurs. At nadir, the illuminated area is expressed by

tic =  é ê ë p 2 + g ù ú û r2 ,,

(76)

where g = sin^-1 [[(tan(a))/(tan(q_s))]] [19]. This is valid for all azimuth angles f. On the backscattering side of the principal plane (Fig. 19a), the shadow can be seen while q_v < q_s. For q_v £a we have

tic = yE ·(xE-xB) + rxA ·[ sin(g)cos(g) -g+p/2 ]

(77)

where (x_E,y_E) denotes point E. x_E = x_A ·sin(g+f) and y_E = r·cos(g+f) As q_v increases, the triangular part appears to the viewer. Two shadowed areas can be seen on the triangle part. Outside the principal plane, the symmetry does not exist, but the same equations can still be used. The third plate of Fig. 19b shows such a case. The computation of the shaded area is an integration of areas inside the ellipse which is separated into four parts to facilitate the selection of integral limits. Figure 19b for a > q_v > 0 we have the shadowed area in three parts of the ellipse, it is calculated as follow: 

ó õ yE 0 é ê ë y-bBE mBE - xa r   ______ Ör2- y2 ù ú û dy  +  p 4 ·r2 cos(qv) - |xB·bBF| 2 + ó õ xF 0 é ê ë mBFx+bBF- r xa   _______ ÖxA2- x2 ù ú û dy.

(78)

The first term is for the area formed by the points E, B, and A, denoted EBA, the second and third are for the area ABGC and GFC, respectively. The point G is the intercept of [`BF] on y. The area illuminated (t_ic) is found by substracting ( 78) from (75).

In (78) m_BE and b_BE are the slope and intercept of [`BE] on y defined by the equation of the curve that has the 2 points B and E. m_BF and b_BF have similar definitions. y_E and y_F do not depend on q_v. Once they are found at nadir, their values can be apply to other view angle at any given f. When q_v > a, the amount of shadowed surface in the triangle is calculated as follows: the integration is like the one in (75). We divide the triangle into two parts separated by the x-axis. The areas of the two triangles are:

ó õ yI 0 é ê ë xA r   ______ Ör2- y2 - y-bBE mBE ù ú û dy

(79)

and

ó õ yJ 0 é ê ë xA r   ______ Ör2- y2 - y-bBF mBF ù ú û dy.

(80)

The shadow on the ellipse is, as above, calculated by separating the ellipse into four parts: HOAJ, GOAI, FGOK, and EHOK.

more text

this is more text