Modeling the Directional Clumping Index of Crop and Forest

The Clumping Index (Ω) was introduced to quantify the spatial distribution pattern of vegetation elements. It is crucial to improve the estimation accuracy of vital vegetation parameters, such as Leaf Area Index (LAI) and Gross Primary Production (GPP). Meanwhile, the parameterization of Ω is challenging partly due to the varying observations of canopy gaps from different view angles. Many previous studies have shown the increase of Ω with view zenith angle through samples of gap size distribution from in situ measurements. In contrast, remote sensing retrieval algorithms only assign a constant value for each biome type to roughly correct the clumping effect as a compromise between the accuracy and efficiency. In this paper, analytical models are proposed that estimate the directional clumping index (Ω(θ)) of crop and forest at canopy level. The angular variation trend and magnitude of crop Ω(θ) was analyzed within row structure where vegetation elements are randomly spaced along rows. The forest model predicts Ω(θ) with tree density, distribution pattern, crown shape, trunk size, and leaf area and angle distribution function. The models take into account the main directional characteristics of clumping index using easy-to-measure parameters. Test cases showed that Ω(θ) magnitude variation for black spruce forest was 102.3% of the hemispherical average clumping index (Ω̃), whereas the Larch forest had 48.7% variation, and row crop variation reached 32.4%. This study provided tools to assess Ω(θ) of discontinuous canopies.


Introduction
Canopy clumping affects the gap fraction and redistributes radiation interception within the canopy [1][2][3][4].The assumption of random distributions causes underestimation of leaf area index (LAI) as compared to considering clumping, since the clumping phenomenon of natural forest increases the gap fraction [5][6][7].Therefore, the effective LAI that was indirectly obtained from gap fraction measurements needs to be interpreted together with the clumping condition to approach the actual LAI, i.e., the one-sided leaf area per unit horizontal ground surface area [8][9][10][11].
Due to this complexity, the study and application of Ω faces accuracy-efficiency-robustness trade-off [19].In field campaigns, Ω is derived from gap size statistics or gap fraction, as measured by many ground based optical instruments, such as the Tracing Radiation and Architecture of Canopies (TRAC) system [18] or fish-eye photographs [24].In situ measured Ω can accurately depict the site-specific Ω variation, but lack efficiency.Nevertheless, it is risky to extrapolate the values or relationships beyond the site.
In operational retrieval of MODIS LAI product, the algorithm assigns a single clumping factor to each biome type [25].This is the most simplified method but it guarantees the robustness.Chen et al. [4] developed a satellite level regression model for hemispherical average clumping index from an angular index that is called the normalized difference between hotspot and dark spot (NDHD), which made it possible to retrieve global time series average clumping index [9,26,27].The NDHD method balances the need for computational efficiency while including the pixel specific information.However, the directional variation of clumping index has not been included.
The angular variation of directional clumping index (Ω(θ)) contains abundant vegetation structure information [28].Discontinuous crop and forest are two representative examples.The directional variation pattern of the clumping effect reflects various vegetation canopy morphologies [29,30].
Furthermore, the directionality of clumping index is necessary to be taken into account in vegetation monitoring for two reasons.One reason is the observed clumping varies along with the view direction, which results in the uncertainty in using single clumping index value.Another reason is the accuracy of Ω estimation has an apparent influence on LAI and Gross Primary Production (GPP).Ω estimation error directly and linearly affects LAI retrievals.Normally, LAI measurements from optical instruments and sensors are referred to as effective LAI (LAI e ), the product of Ω and actual LAI (LAI) (half the total leaf area) (LAI e = Ω•LAI).Garrigues et al. [31] tested LAI retrieval sensitivity while using digital hemispherical photography over 10 crop types.The clumping index induced difference in LAI estimation at 30-m scale was as high as 31% (relative root-mean-square error or relative RMSE).Weiss [14] suggested that the accurate estimation of clumping index for each zenith angle will increase LAI estimation accuracy.Piayda et al. [32] applied angularly dependent leaf clumping index in determining LAI, which led to a 30% higher estimation when compared to the indirect measurements using LAI-2000, but matched precisely with the direct measurements using litter traps.Chen et al. [33] tested the effects of ignoring foliage clumping characteristics in global GPP mapping, and found that neglecting clumping index will increase the uncertainty in GPP estimation, particularly in tropical regions.Using accurate LAI, but ignoring clumping index, causes an overestimation of global GPP by 12%; applying the effective LAI without considering the clumping effect would underestimate global GPP by 9%.Comparing the effect of clumping index on LAI/GPP estimation and the directional variation magnitude of Ω [2,21], it is necessary to apprehend that if angular variation is not provided as well as the average, large errors can occur in the inverted leaf angle distribution and LAI for canopies with direction dependent clumping [34].
Many studies have discussed the directional variation of Ω using the gap size distribution measurements from in situ instruments, such as TRAC or digital hemispherical photography.The efficiency in data measurement is the bottleneck in wide application of the gap size distribution model.Kuusk [22] developed a simulation model that had extended the understanding of directional clumping effect, but it was difficult to be applied for practical calculations, as many input parameters cannot be measured.Kucharik [21] modeled the directional clumping index using an empirical equation to extend its availability for more species, although the accuracy of coefficients is restricted in the regression hull.This paper aims to propose analytical models to estimate directional clumping index of crop and forest at canopy level.
According to the Beer's law [35], the clumping index can be estimated as a function of gap fraction of vegetation canopy.Crop canopy illustrates an apparent clumping index angular dependence due to its gap fraction variation in discrete periodic row structures [22,36,37].Fan et al. [37] proposed an analytical model for the crop clumping index using the directional average gap fraction, where each hedgerow was regarded as the integration of four parts with different average LAI.This assumption led to an accurate expression of crop clumping index, although the integration component is incorrect in concept.Peng et al. [38] rectified the method and applied the directional clumping index in retrieving crop albedo.The average gap fraction of forest was estimated according to the interaction between light regime and the canopy, such as the widely used Geometric-Optical (GO) model that was proposed by Li et al. [39], Nilson [7], and Xu [40].In the GO model, forest gap fraction is the complement of the projected area fraction on the ground.Nilson [7] proposed the gap probability was determined by the number of trunks occurring in the shadowed region of a single tree.It is assumed each tree crown was ellipsoid, which fits most practical cases.The measurements and simulations [21] suggest that the largest gaps between crowns are a major cause of the canopy clumping effect.
This paper proposed analytical models for crop and forest directional clumping index based on gap fraction estimation.To extend the applicability of the model, we emphasized the angular variation of Ω due to the first-level structure, i.e., the row structure of crop, and the abundance of crown shapes coupled with complex tree spatial distribution of forest.
Section 1 introduced the importance of clumping index and the necessity of directional clumping index according to the previous studies.Section 2 introduces the derived analytical model of angular clumping index variation for crop and forest.Section 3 demonstrates validation of the model, while using in situ results to validate the reliability of the proposed methods, and we evaluated the variation range of Ω(θ).Section 4 analyzes the influence of various parameters on directional clumping index using the proposed methods.Section 5 summarizes the main conclusion and proposes the limitations in the model and analysis.The symbols used in this paper are listed in Table 1.

Modeling and Evaluation
Crop Ω generally varies with azimuth angle, whereas forest Ω is generally symmetric in the azimuth dimension.Different algebraic models are firstly proposed for their respective directional clumping index estimation.
Simulation is a practical choice to validate algebraic expression.The Monte Carlo (MC) method was employed for crops [38].The morphological statistics was proposed for forests, which quantifies the gap fraction through automatic counting that is based on virtual forest scenarios.Simulated results can capture the general variation trends of directional gap fractions and provide reference results.

Crop Model
We regard the crop as clumps of vegetation elements at the canopy level [37].The vegetation elements are randomly distributed within hedgerows [41].Following the Beer's law, the average gap fraction of the vegetation canopy is a function of view zenith angle and its structural parameters, including clumping index, LAI, and leaf angle distribution [42], where p(θ, ϕ) is the gap fraction; Ω(θ) refers to the directional clumping index; G θ is a parameter to describe canopy structure; µ is the cosine of zenith angle θ; and, LAI a denotes the spatially-averaged leaf area index.
The average gap fraction of hedgerow crop canopy can be estimated using Yan et al.'s algorithm [41].
where A is a row unit consisting of canopy width (A 1 = a 1 / sin ϕ) and row spacing (A 2 = a 2 / sin ϕ) at a relative azimuth angle of ϕ with the row direction; a 1 and a 2 represents the canopy width and row spacing that vertical to the row direction.H is the height of the crop.Then, Ω(θ) can be calculated from the average gap fraction, p(θ, ϕ), i.e., where the equivalent LAI (LAI a ), includes the influence of foliage and non-foliage elements.Kucharik et al. [16,21] proposed that the canopy clumping index is a function of path length through the canopy.More specifically, it is a function of the effective path length in proportion to the encountered vegetation elements, as it is related to the energy attenuation in the canopy.For crop canopy, this effective path length does not count the row gap.Therefore, the maximum effective path length shows ups and downs while penetrating more rows and row gaps at increasing zenith angles, which is illustrated as the line A-E in Figure 1.
Accordingly, the variation of Ω(θ) with zenith angle can be pinpointed while using the local feature angles, θ ci , corresponding to the peaks of the Ω(θ) curve with zenith angle (Figure 2).The feature angles are determined by the row structure, where i = 1, 2, 3, . . .i denotes the number of rows that photons penetrated through in the crop canopy; and A = A 1 + A 2 = a 1 / sin ϕ + a 2 / sin ϕ.The subscript c and d corresponds to the local peak and valley value.
Ω(θ) increases along with zenith angle before reaching the first feature angle (Figure 2).In practice, if the fluctuation (normally < 0.05) can be tolerated for Ω(θ), we can regard Ω(θ) as approximately stable after the first turning point θ c .Therefore, θ c can be used to predict the main Ω(θ) trend, which is regarded as the range angle of Ω(θ), feature angles are determined by the row structure, where i = 1, 2, 3, … i denotes the number of rows that photons penetrated through in the crop canopy; and Ω(θ) increases along with zenith angle before reaching the first feature angle (Figure 2).In practice, if the fluctuation (normally < 0.05) can be tolerated for Ω(θ), we can regard Ω(θ) as approximately stable after the first turning point c θ .Therefore, c θ can be used to predict the main Ω(θ) trend, which is regarded as the range angle of Ω(θ), where i = 1, 2, 3, … i denotes the number of rows that photons penetrated through in the crop canopy; and Ω(θ) increases along with zenith angle before reaching the first feature angle (Figure 2).In practice, if the fluctuation (normally < 0.05) can be tolerated for Ω(θ), we can regard Ω(θ) as approximately stable after the first turning point c θ .Therefore, c θ can be used to predict the main Ω(θ) trend, which is regarded as the range angle of Ω(θ),

Forest Model
The derivation of forest Ω(θ) uses the average gap fraction (3), which is referred to Nilson's model [7].We modelled the tree crown as vertical ellipsoid for three reasons.First, the vertical ellipsoid is axially symmetrical, which accords with our assumption that the variation of clumping index along with azimuth angle is not considered.Second, the crowns of different species can be discriminated though the ratio between crown depth and crown width, which is related to the variation of the clumping index along with zenith angle.Third, the ellipsoid can be easily parameterized using crown depth and crown width in the model.
For randomly distributed forest, the gap fraction at canopy level can be expressed as where p i (θ) is the intra-crown directional gap fraction of the tree, h/2R is a crown shape factor that discriminates different trees types or species, and nπR 2 is the crown closure where the overlapped regions of crown projections have been multiply counted, i.e., C CR = nπR 2 .Combining (3) and (4), However, ( 6) and (7) do not consider the trunk, which can affect optical measurements for some high trees with tall or wide trunks [43].The average gap fraction at canopy level, incorporating trunk effects, can be expressed as and the directional Ω(θ) can be expressed as Alternatively, if we count the overlapped regions of crown projections only once in the canopy closure calculation, then C CA = 1 − e −n πR 2 , since it is equivalent to the canopy gap fraction at the nadir direction.Then, the canopy average gap fraction and Ω(θ) can be expressed as, respectively, and if considering the tree trunk, and Equations ( 7) and (9) are Ω(θ) expressions for multiple counted canopy closure, and Equations ( 11) and (13) are Ω(θ) expressions for once counted canopy closure.

Simulation Method of Directional Clumping Index
To evaluate the performance of the analytical model, we used the Monte Carlo simulation methods to obtain cross validation information when considering the ground measurements are limited.Compared to the analytical methods, the Monte Carlo simulations are based on the repeated random sampling and statistical analysis.They provide numerical solutions to the unknown parameters with the same simplifications and assumptions as the analytical methods so that the results can be inter-compared.Here, comparing with current clumping index models is not applicable because of the definition difference.The models aim to obtain the average clumping index Ω rather than directional clumping index Ω(θ).
For crop, we used a Monte Carlo (MC) method to simulate the Ω(θ), as in Section 2.2.1.For forest, we used a virtual forest scenario to simulate the forest gap fraction and Ω(θ) to provide validations for the algebraic model.The detailed setting of the forest scenario is described in Section 2.2.2.The procedure of obtaining canopy-level directional gap fraction is introduced in Section 2.2.3.This simulation mainly reflects the inter-crown gap fraction variation, as this dominantly causes the gap size distribution function to deviate from random at canopy level [19,21].

Crop Directional Clumping Index Simulation
The MC method simulated the radiative transfer process of photons in the canopy by generating the incident photons, determining the scenarios regarding photons colliding with leaves or soil, randomly assigning the scattering direction, calculating the free paths of the photons, and counting the photons satisfying specific conditions.The detailed procedures and parameter setting was introduced in [38].
A set of parameters from maize canopy was set as the hypothesis example (a 1 = 0.2, a 2 = 0.3, H = 1, r h g = 0.22, ρ = 0.0999, T = 0.1303, Table 1).We calculated Ω(θ, ϕ) while using the MC method for the same scenario.In the MC method, the radiative photon transfer process in the canopy is simulated for predefined canopy structural parameters.Each photon was traced and collisions were recorded.The statistics from many photon behaviors approximately reflect actual canopy characteristics and provide validation.Therefore, the gap fraction is the fraction of photons penetrating the canopy, and Ω(θ, ϕ) can be calculated from Equation (3).

Forest Scenario Settings
Our goal was to directly count the gap fraction by traversing, so we must first construct a forest scenario.The main problems involve how to depict each tree and how to define the distribution pattern for multiple trees.The scenes can be extended to be infinite through juxtaposing identical scenes.
The basic tree crown shape is modelled as a generalized ellipsoid as in the analytical forest model (Section 2.1.2).Then the crown shape could be defined by the major and minor axis lengths.The trunk was treated as a cylinder defined by height and radius at breast height.We defined similar trees in each simulation for efficiency, but this could be easily extended to various tree shapes by introducing random offsets.Then, we define the tree positions according to the different distribution patterns.If the trees are distributed randomly, the distribution follows Poisson's law [8,44,45].
However, trees are likely to cluster in different sized groups rather than truly distribute randomly.Chen and Leblanc [8] assumed the distribution of the groups and the trees in each group both followed Poisson's law.This dual Poisson or Neyman distribution fits the actual situation better than the simple Poisson's law [45][46][47].The Neyman distribution PDF can be expressed as where m 2 = 1 implies that the tree distribution is close to, but not the same as, the random case.
Canopy closure is directly linked to the gap fraction.Thus, the simulation was conducted using increasing m (Poisson distribution) and m 1 (Neyman distribution), corresponding to increasing stand density (trees/m 2 ).

Forest Directional Clumping Index Simulation
The whole forest canopy was segmented into small inclined cylinders along the view direction, while considering the leaf distribution heterogeneity (Figure 3).Average LAI can be calculated for each cylinder according to the tree distribution pattern, i.e., the predefined leaf volume density and cylinder height (related to view zenith angle, θ).Assuming leaves distribute homogeneously within each cylinder, Ω(θ) is 1 for each cylinder at any zenith angle.The cylinder gap fraction p i can be calculated according to Beer's law as a function of u, where the canopy height is converted to cylinder height through cos θ, and G θ = 0.5.The average gap fraction, p(θ), of the whole canopy can be calculated as the average over all cylinders, and Ω can be obtained from (3).

Forest Directional Clumping Index Simulation
The whole forest canopy was segmented into small inclined cylinders along the view direction, while considering the leaf distribution heterogeneity (Figure 3).Average LAI can be calculated for each cylinder according to the tree distribution pattern, i.e., the predefined leaf volume density and cylinder height (related to view zenith angle, θ).Assuming leaves distribute homogeneously within each cylinder, Ω(θ) is 1 for each cylinder at any zenith angle.The cylinder gap fraction i p can be calculated according to Beer's law as a function of u, where the canopy height is converted to cylinder height through cosθ , and 0.5 . The average gap fraction, ( ) p θ , of the whole canopy can be calculated as the average over all cylinders, and Ω can be obtained from (3).Key parameters for the simulation include the average LAI of each tree, tree height, crown depth, crown radius, tree density, and parameters defining the tree distribution pattern and research area size.

Forest Ground Measurements
In this paper, two groups of forest Ω(θ) ground measurements were used to validate the analytical model.The ground measurements were collected from published field campaigns and new in situ experiments.Simulation parameters for forest characteristics were set according to  Key parameters for the simulation include the average LAI of each tree, tree height, crown depth, crown radius, tree density, and parameters defining the tree distribution pattern and research area size.

Forest Ground Measurements
In this paper, two groups of forest Ω(θ) ground measurements were used to validate the analytical model.The ground measurements were collected from published field campaigns and new in situ experiments.Simulation parameters for forest characteristics were set according to practical forest features.The first group of data was obtained from an old Black Spruce (Picea mariana) forest, near Candle Lake, Saskatchewan [8,18] while using TRAC with a measurement error of about 5%.
A second group of measurements was collected at the Saihanba experimental site within Hebei Province in China (117.32 • E, 42.40 • N).The domain was covered by steppe and forest steppe, with Larch (Larix gmelinii (Rupr.)Rupr.), one of the main tree types.
We measured the Larch gap fraction using hemispherical photographs obtained from a Canon EOS 50D camera with a fish-eye lens, providing orthographic projection and a 180 • field of view (Figure 4).The camera setting is introduced in [48].53 images were taken within the 45 m × 45 m square.The upward looking hemispherical photograph was further processed using Hemiview 2.1 software (http://www.delta-t.co.uk/product/hemiview/, 5 June 2017) to calculate canopy multiple angle gap fraction.Each hemispherical image was divided into 18 concentric rings between 0 • and 90 • zenith angle, with a zenith resolution of 5 • .The directional gap fraction is calculated for each zenithal ring.The measurement error mainly depends on the threshold setting to discriminate vegetation from sky conditions, which is about 1-3%.Simulation parameters are listed in Table 2, column 2, as collected from the field measurements.
A second group of measurements was collected at the Saihanba experimental site within Hebei Province in China (117.32°E,42.40°N).The domain was covered by steppe and forest steppe, with Larch (Larix gmelinii (Rupr.)Rupr.), one of the main tree types.
We measured the Larch gap fraction using hemispherical photographs obtained from a Canon EOS 50D camera with a fish-eye lens, providing orthographic projection and a 180° field of view (Figure 4).The camera setting is introduced in [48].53 images were taken within the 45 m × 45 m square.The upward looking hemispherical photograph was further processed using Hemiview 2.1 software (http://www.delta-t.co.uk/product/hemiview/, 5 June 2017) to calculate canopy multiple angle gap fraction.Each hemispherical image was divided into 18 concentric rings between 0° and 90° zenith angle, with a zenith resolution of 5°.The directional gap fraction is calculated for each zenithal ring.The measurement error mainly depends on the threshold setting to discriminate vegetation from sky conditions, which is about 1-3%.Simulation parameters are listed in Table 2, column 2, as collected from the field measurements.

Crop Directional Clumping Index Validation
The variation of Ω(θ, ϕ) with View Azimuth Angle (ϕ) under different θ is illustrated in Figure 5.The differences between the modeled and simulated results (Figure 6) were randomly distributed with an average of 0, and maximum relative discrepancy <1.5%.The MC simulated and algebraic Ω(θ, ϕ) are strongly consistent, which provides confidence to analyze crop Ω(θ, ϕ) variation while using the algebraic model.
decreases, thus c θ grows; when increases from 90° to 150°, A increases, then c θ drops.When the is 0° or 180°, i.e., along with the row direction, the gap fraction of row plant is obviously larger than other directions as the row gap always occupies a fraction of the field of view.The Ω( , ) decreases with the growing view zenith angle, which means the crop departures more away from the random status.On the contrary, when the azimuth angle is larger, the Ω( , ) starts growing with view zenith angle and reaches a stable value around one after the view zenith angle higher than the range angle.When the azimuth angle reaches 90°, the ratio of equivalent row width to row height gets its minimum value so that the range angle is smaller than other azimuth directions.These variation trends verify the patterns of crop Ω( , ) in Section 2.1.1 and they show the crop Ω( , ) dispersion at different view angles from constant status.Other simulation parameters: decreases, thus c θ grows; when increases from 90° to 150°, A increases, then c θ drops.When the is 0° or 180°, i.e., along with the row direction, the gap fraction of row plant is obviously larger than other directions as the row gap always occupies a fraction of the field of view.The Ω( , ) decreases with the growing view zenith angle, which means the crop departures more away from the random status.On the contrary, when the azimuth angle is larger, the Ω( , ) starts growing with view zenith angle and reaches a stable value around one after the view zenith angle higher than the range angle.When the azimuth angle reaches 90°, the ratio of equivalent row width to row height gets its minimum value so that the range angle is smaller than other azimuth directions.These variation trends verify the patterns of crop Ω( , ) in Section 2.1.1 and they show the crop Ω( , ) dispersion at different view angles from constant status.The different shapes of Ω(θ, ϕ) − ϕ curves reflect the effect of the range angle θ c of Ω(θ, ϕ) pattern.When ϕ increases from 30 • to 90 • (perpendicular to the row direction), row period width A decreases, thus θ c grows; when ϕ increases from 90 • to 150 • , A increases, then θ c drops.When the ϕ is 0 • or 180 • , i.e., along with the row direction, the gap fraction of row plant is obviously larger than other directions as the row gap always occupies a fraction of the field of view.The Ω(θ, ϕ) decreases with the growing view zenith angle, which means the crop departures more away from the random status.On the contrary, when the azimuth angle is larger, the Ω(θ, ϕ) starts growing with view zenith angle and reaches a stable value around one after the view zenith angle higher than the range angle.When the azimuth angle reaches 90 • , the ratio of equivalent row width to row height gets its minimum value so that the range angle is smaller than other azimuth directions.These variation trends verify the patterns of crop Ω(θ, ϕ) in Section 2.1.1 and they show the crop Ω(θ, ϕ) dispersion at different view angles from constant status.

Forest Directional Clumping Index Validation
Figure 7 shows the comparison result in the old Black Spruce forest.Since the simulated and measured Ω(θ) were both calculated from the gap fraction, we directly compared the gap fraction for validation.There is high consistency between the measured and simulated gap fraction, validating the proposed simulation method.At zenith angles of >50 • , the simulated gap fraction tends to be slightly lower than the measured value due to multiple photon scattering impacting on optical measurements [8,18].

Forest Directional Clumping Index Validation
Figure 7 shows the comparison result in the old Black Spruce forest.Since the simulated and measured Ω(θ) were both calculated from the gap fraction, we directly compared the gap fraction for validation.There is high consistency between the measured and simulated gap fraction, validating the proposed simulation method.At zenith angles of >50°, the simulated gap fraction tends to be slightly lower than the measured value due to multiple photon scattering impacting on optical measurements [8,18].2.
We calculated the gap fraction using Equation (12), with the intra-crown gap pre-processed through simulation as a Look Up Table .The other parameters, including tree density, crown depth, and radius, were set as for previous simulation scenarios.The algebraic method includes the ideal assumption that the trees were Neyman distributed.However, the simulated Poisson distribution may demonstrate higher clustering than ideal situations [49], due to restrictions of the pseudo random number generator in the simulation.One subset of our simulated tree distribution image is demonstrated in Figure 8.Therefore, the calculated gap fraction from the simulation is a little lower than the expected value.2.
We calculated the gap fraction using Equation ( 12), with the intra-crown gap pre-processed through simulation as a Look Up Table .The other parameters, including tree density, crown depth, and radius, were set as for previous simulation scenarios.The algebraic method includes the ideal assumption that the trees were Neyman distributed.However, the simulated Poisson distribution may demonstrate higher clustering than ideal situations [49], due to restrictions of the pseudo random number generator in the simulation.One subset of our simulated tree distribution image is demonstrated in Figure 8.Therefore, the calculated gap fraction from the simulation is a little lower than the expected value.

Forest Directional Clumping Index Validation
Figure 7 shows the comparison result in the old Black Spruce forest.Since the simulated and measured Ω(θ) were both calculated from the gap fraction, we directly compared the gap fraction for validation.There is high consistency between the measured and simulated gap fraction, validating the proposed simulation method.At zenith angles of >50°, the simulated gap fraction tends to be slightly lower than the measured value due to multiple photon scattering impacting on optical measurements [8,18].2.
We calculated the gap fraction using Equation ( 12), with the intra-crown gap pre-processed through simulation as a Look Up Table .The other parameters, including tree density, crown depth, and radius, were set as for previous simulation scenarios.The algebraic method includes the ideal assumption that the trees were Neyman distributed.However, the simulated Poisson distribution may demonstrate higher clustering than ideal situations [49], due to restrictions of the pseudo random number generator in the simulation.One subset of our simulated tree distribution image is demonstrated in Figure 8.Therefore, the calculated gap fraction from the simulation is a little lower than the expected value.Figure 9 shows the comparison result using the measurements of Larch stand.It shows high consistency between the simulated and measured gap fractions, with small fluctuations at small angles due to random errors and sampling restrictions.When zenith angle of >70 • , discrepancies mainly arise due to the hemispherical image being influenced by topography at large angles.Overlapping between trees in the simulation and approximation of the projection area to accumulated discrete square units determine that the average gap fraction in the simulation method is slightly larger than the model in nadir direction, although the number of trees is equal.This difference weakens with increased zenith angle and disappears between 40 • -60 • , then increases to 0.02 at 80 • .Figure 8. Simulated distribution pattern of black spruce (Neyman distribution, m2 = 4, 1000 trees distributed in 4 ha).The white circles represent the nadir projection of the crowns to mark the tree positions, while the black represents the soil background.
Figure 9 shows the comparison result using the measurements of Larch stand.It shows high consistency between the simulated and measured gap fractions, with small fluctuations at small angles due to random errors and sampling restrictions.When zenith angle of >70°, discrepancies mainly arise due to the hemispherical image being influenced by topography at large angles.Overlapping between trees in the simulation and approximation of the projection area to accumulated discrete square units determine that the average gap fraction in the simulation method is slightly larger than the model in nadir direction, although the number of trees is equal.This difference weakens with increased zenith angle and disappears between 40°-60°, then increases to 0.02 at 80°.Ω(θ) was indirectly validated through the measured, simulated, and modeled gap fraction, and demonstrated that the model is able to predict the main variation characteristics of Ω(θ).  2.

Assessing the Variation Magnitude of Ω(θ)
Using the proposed models, we assessed the relative variation magnitude of Ω(θ) to indirectly reflect its potential influence on LAI/GPP estimation.
We selected three scenarios, one representing a typical crop example, and the other two matching forest measurements.We first calculated the hemispherical average clumping index (Ω ) for each group by angular integration, where Ω represents the commonly used single value of clumping index.The uncertainty due to using Ω rather than angle specific Ω(θ) can be assessed from the position of Ω in the range of Ω( ) (Figure 10).For crop canopies, the discrepancy between Ω( ) and Ω is more obvious before the zenith angle reaches the range angle.Since the range angle of Ω( )~ curve is smaller than 20°, the marker of Ω is out of the box.For forest canopies, the angular variation of the black spruce forest is more intense than the larch forest, due to the different Ω(θ) variation range determined by their tree species and growth status.However, all of the variation amplitude relative to Ω are large and they should not be neglected.
To further quantify the potential uncertainty caused by ignoring variation of Ω( ) with zenith angle, we estimated the relative angular variation magnitude of Ω( ) using the index where AV values for the example clumping index datasets are listed in Table 3. Ω( ) variation amplitude for crops can reach 32.4% of its average value.Black spruce has tall and slender crowns, i.e., larger crown shape factor (Section 4.2.4), which leads to larger Ω( ) variation range, with Av reaching 102.3%.Larch forest canopy AV is somewhat lower, mainly due to its relatively smaller Ω( ) variation range related to its crown shape factor.This demonstrates that potential uncertainty  2.
Ω(θ) was indirectly validated through the measured, simulated, and modeled gap fraction, and demonstrated that the model is able to predict the main variation characteristics of Ω(θ).

Assessing the Variation Magnitude of Ω(θ)
Using the proposed models, we assessed the relative variation magnitude of Ω(θ) to indirectly reflect its potential influence on LAI/GPP estimation.
We selected three scenarios, one representing a typical crop example, and the other two matching forest measurements.We first calculated the hemispherical average clumping index ( Ω) for each group by angular integration, where Ω represents the commonly used single value of clumping index.The uncertainty due to using Ω rather than angle specific Ω(θ) can be assessed from the position of Ω in the range of Ω(θ) (Figure 10).For crop canopies, the discrepancy between Ω(θ) and Ω is more obvious before the zenith angle reaches the range angle.Since the range angle of Ω(θ) ∼ θ curve is smaller than 20 • , the marker of Ω is out of the box.For forest canopies, the angular variation of the black spruce forest is more intense than the larch forest, due to the different Ω(θ) variation range determined by their tree species and growth status.However, all of the variation amplitude relative to Ω are large and they should not be neglected.
To further quantify the potential uncertainty caused by ignoring variation of Ω(θ) with zenith angle, we estimated the relative angular variation magnitude of Ω(θ) using the index where A V values for the example clumping index datasets are listed in Table 3. Ω(θ) variation amplitude for crops can reach 32.4% of its average value.Black spruce has tall and slender crowns, i.e., larger crown shape factor (Section 4.2.4), which leads to larger Ω(θ) variation range, with Av reaching 102.3%.Larch forest canopy A V is somewhat lower, mainly due to its relatively smaller Ω(θ) variation range related to its crown shape factor.This demonstrates that potential uncertainty can arise by ignoring Ω(θ) variation in practice.Comparing the relative variation magnitude with the allowable error range will help to determine whether or not Ω(θ) should be treated as constant in practical situations.
Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 22 can arise by ignoring Ω( )variation in practice.Comparing the relative variation magnitude with the allowable error range will help to determine whether or not Ω( ) should be treated as constant in practical situations.
Another main application for Ω( ) is to estimate the quantity of the terrestrial carbon sink and its dynamic variability.One representative parameter is GPP [33,50].The proposed method can be used in GPP estimation to explore whether it would improve the accuracy, which would be a fruitful future research direction.

Sensitivity Analysis
According to the analytical models, the factors influencing crop Ω(θ) include LAI, leaf foliage distribution, and hedgerow structure, whereas tree density, tree distribution, crown shape, and trunk size influence the forest clumping index.The contribution of these factors are analyzed in this section to assess their relative importance.

Leaf Area Index
Figure 11 shows that Ω(θ) decreases with increasing LAI in all directions, since heterogeneity is more apparent for larger intra-ridge LAI.Before the zenith angle reaches the range angle, Ω(θ) varies greatly with crop canopy LAI.The view zenith angle of satellite observations is mainly smaller than the range angle, e.g., MODIS view zenith angle is between 0° and 65°.Therefore, Ω(θ) is essential to LAI estimation while using satellite data.Another main application for Ω(θ) is to estimate the quantity of the terrestrial carbon sink and its dynamic variability.One representative parameter is GPP [33,50].The proposed method can be used in GPP estimation to explore whether it would improve the accuracy, which would be a fruitful future research direction.

Sensitivity Analysis
According to the analytical models, the factors influencing crop Ω(θ) include LAI, leaf foliage distribution, and hedgerow structure, whereas tree density, tree distribution, crown shape, and trunk size influence the forest clumping index.The contribution of these factors are analyzed in this section to assess their relative importance.

Leaf Area Index
Figure 11 shows that Ω(θ) decreases with increasing LAI in all directions, since heterogeneity is more apparent for larger intra-ridge LAI.Before the zenith angle reaches the range angle, Ω(θ) varies greatly with crop canopy LAI.The view zenith angle of satellite observations is mainly smaller than the range angle, e.g., MODIS view zenith angle is between 0 • and 65 • .Therefore, Ω(θ) is essential to LAI estimation while using satellite data.

Leaf Angle Distribution
Leaf angle distribution (G θ ) is a key parameter influencing Ω(θ) estimation.We estimated crop Ω(θ) for different G θ functions (Erectophile, Plagiophile, Spherical, Extremophile, and Planophile) (Figure 12) [51].Figure 13 shows that the effect of different leaf angle distribution is more evident closer to nadir direction.For zenith angle beyond the θ c , this effect lessens as the crop canopy approaches continuous.Leaf angle distribution and Ω(θ) influence accurate LAI estimation from satellite observations with smaller zenith angles than the range angle.Leaf angle distribution ( G θ ) is a key parameter influencing Ω(θ) estimation.We estimated crop Ω(θ) for different G θ functions (Erectophile, Plagiophile, Spherical, Extremophile, and Planophile) (Figure 12) [51].Figure 13 shows that the effect of different leaf angle distribution is more evident closer to nadir direction.For zenith angle beyond the θc, this effect lessens as the crop canopy approaches continuous.Leaf angle distribution and Ω(θ) influence accurate LAI estimation from satellite observations with smaller zenith angles than the range angle.Row structure is characterized by canopy height, H, and row period width, A. As discussed above, Ω(θ) increases with zenith angle until the range angle, which is controlled by the row structure.Figure 14 shows calculated Ω(θ) for various row structures that are discriminated by A/H, where H = 1, and 1 a / 2 a = 2:3.Thus, when A/H increases with growing A, the range angle grows, which is consistent with (5).Generally, the LAI and A/H both have apparent influence on Ω(θ), while has less influence than the other two factors.However, LAI and A/H are inter-correlative, as we assume that the leaf volume density keeps constant [41], and their influence will partly compensate each other.The influence of all three factors are less obvious at large zenith angles, because Ω(θ) reaches a relatively stable value close to 1 after the zenith angle exceeds the range angle, when the discontinuous crop canopy can be approximately regarded as continuous vegetation canopy.

Row Structure
Row structure is characterized by canopy height, H, and row period width, A. As discussed above, Ω(θ) increases with zenith angle until the range angle, which is controlled by the row structure.Figure 14 shows calculated Ω(θ) for various row structures that are discriminated by A/H, where H = 1, and a 1 /a 2 = 2:3.Thus, when A/H increases with growing A, the range angle grows, which is consistent with (5).Row structure is characterized by canopy height, H, and row period width, A. As discussed above, Ω(θ) increases with zenith angle until the range angle, which is controlled by the row structure.Figure 14 shows calculated Ω(θ) for various row structures that are discriminated by A/H, where H = 1, and 1 a / 2 a = 2:3.Thus, when A/H increases with growing A, the range angle grows, which is consistent with (5).Generally, the LAI and A/H both have apparent influence on Ω(θ), while has less influence than the other two factors.However, LAI and A/H are inter-correlative, as we assume that the leaf volume density keeps constant [41], and their influence will partly compensate each other.The influence of all three factors are less obvious at large zenith angles, because Ω(θ) reaches a relatively stable value close to 1 after the zenith angle exceeds the range angle, when the discontinuous crop canopy can be approximately regarded as continuous vegetation canopy.Generally, the LAI and A/H both have apparent influence on Ω(θ), while G θ has less influence than the other two factors.However, LAI and A/H are inter-correlative, as we assume that the leaf volume density keeps constant [41], and their influence will partly compensate each other.The influence of all three factors are less obvious at large zenith angles, because Ω(θ) reaches a relatively stable value close to 1 after the zenith angle exceeds the range angle, when the discontinuous crop canopy can be approximately regarded as continuous vegetation canopy.

Factors Influencing Forest Clumping Index Variation
Crown coverage and crown depth dominates the directional gap fraction and Ω(θ) [21].In natural cases, crown coverage not only relies on tree density, but also the distribution patterns.Morphological characteristics of tree species also influence the gap fraction.Therefore, we investigated the influence of tree density, tree distribution, crown shape, and trunk size on the clumping index.When we design the simulation experiment and select the parameter values, the following factors were considered.
(1) Comparability.The results among different parameters' results are generally comparable.Some basic setting of the scenarios keeps the same as the ground measured parameters in the black spruce forest; and, (2) Representativeness.For example, we choose an average crown shape of deciduous forest that the crown depth is three times of the crown radius [52].To avoid the violation to the comparability, the total tree volume was kept constant, so the crown depth was set to 2.28 m and the crown radius was 0.76 m.

Tree Density
Tree density was varied with individual tree parameters held constant (Table 2, Column 3), and Ω(θ) estimated for each case, as shown in Figure 15.Forest canopy Ω(θ) increases with growing view zenith angle, but since the distribution of vegetation elements is more heterogeneous in forest than crop canopies, Ω(θ) at large zenith angle cannot reach 1, although it also achieves relative stability.
following factors were considered.(1) Comparability.The results among different parameters' results are generally comparable.Some basic setting of the scenarios keeps the same as the ground measured parameters in the black spruce forest; and, (2) Representativeness.For example, we choose an average crown shape of deciduous forest that the crown depth is three times of the crown radius [52].To avoid the violation to the comparability, the total tree volume was kept constant, so the crown depth was set to 2.28 m and the crown radius was 0.76 m.

Tree Density
Tree density was varied with individual tree parameters held constant (Table 2, Column 3), and Ω(θ) estimated for each case, as shown in Figure 15.Forest canopy Ω(θ) increases with growing view zenith angle, but since the distribution of vegetation elements is more heterogeneous in forest than crop canopies, Ω(θ) at large zenith angle cannot reach 1, although it also achieves relative stability.
For forests with a given crown coverage (in this simulation, approximately > 0.1), inter-crown gaps disappear rapidly with growing zenith angle, as soil in the field of view is gradually replaced by neighboring trees at the oblique angle.Therefore, Ω(θ) increases for 0-50°.In contrast, for sparse forest, e.g., crown coverage = 0.031 in Figure 15, inter-crown gap dominants the field of view and it cannot be totally counteracted through increasing zenith angle.Thus, Ω(θ) increase for sparse forest is not as apparent as for denser forest and the clumping index generally increases with tree density.An obvious difference can be observed between the N = 1000 line and other lines due to the mutual shadowing and overlapping among crowns.Furthermore, Ω(θ) values for a given zenith angle do not change monotonically with tree density, but they show fluctuations due to changes in the tree arrangement for simulation scenarios, since tree positions were randomly determined.This phenomenon also can be seen in Figures 16-18.However, the main trend of Ω(θ) variation remains similar across these factors.For forests with a given crown coverage (in this simulation, approximately > 0.1), inter-crown gaps disappear rapidly with growing zenith angle, as soil in the field of view is gradually replaced by neighboring trees at the oblique angle.Therefore, Ω(θ) increases for 0-50 • .In contrast, for sparse forest, e.g., crown coverage = 0.031 in Figure 15, inter-crown gap dominants the field of view and it cannot be totally counteracted through increasing zenith angle.Thus, Ω(θ) increase for sparse forest is not as apparent as for denser forest and the clumping index generally increases with tree density.An obvious difference can be observed between the N = 1000 line and other lines due to the mutual shadowing and overlapping among crowns.Furthermore, Ω(θ) values for a given zenith angle do not change monotonically with tree density, but they show fluctuations due to changes in the tree arrangement for simulation scenarios, since tree positions were randomly determined.This phenomenon also can be seen in Figures 16-18.However, the main trend of Ω(θ) variation remains similar across these factors.
clumping effect between trees and intra-crown gap fraction increase.However, this combined effect on Ω(θ) is not large.When m2 ≤ 5, Ω(θ) difference between sample lines <0.02.This demonstrates that the slight non-random distributions in natural forests has little influence on canopy level Ω(θ).But, if the non-random distribution is dominant, its influence on Ω(θ) cannot be ignored.For example, the Ω(θ) difference between the sample line with m2 = 5 and that with m2 = 20 can reach 0.08 at large zenith angles.

Trunk
The proposed algebraic models for forest clumping index estimation parameterized tree trunk influence.The influence of tree trunk on gap fraction would be great for high trunk tree species, such as the black spruce forest that is considered above.
Figure 17 shows that the influence of trunk is more important at higher zenith angles, because the fraction of visible trunk grows with zenith angle from only the cross section at nadir direction.Higher trunk fraction in the field of view leads to lower gap fraction, causing increased Ω(θ) for constant average LAI.The variation magnitude depends on the forest scenario, which is larger than 0.1 in this example and it should not be neglected.
Note that increased trunk size would increase the equivalent average LAI, since the trunk is one of the vegetation elements, but the average LAI was set as constant in the simulation (Figure 17) to realize monofactor sensitivity analysis and highlight trunk height influence on directional Ω(θ).

Crown Shape
The tree crown was simplified as an ellipsoid in the proposed model.One convenience of this approach is that different species can be discriminated through the crown shape factor / 2 h R combined with horizontal radius, R.
To highlight crown shape influence on canopy Ω(θ), we kept the crown column constant while changing the crown shape factor (Figure 18).Table 2 shows that crown depth, radius, and their ratio change synchronously.
As the ratio between crown depth and radius increases, there is more obvious increasing Ω(θ) trend with zenith angle.When 2 / 2 h R = , i.e., the crown is spherical, there is only slight Ω(θ) increase with increasing zenith angle (~0.2 in this example).However, when 8 / 2 h R = , implying that a deep canopy, Ω(θ) increased with growing zenith angle, and approaches 1.
At near nadir direction, lower crown depth and radius ratios correspond to larger horizontal projection for each crown when the column is fixed.This larger projection corresponds to lower gap fraction and results in larger Ω(θ).However, curves with higher crown depth and radius ratios demonstrate Ω(θ) increase with growing zenith angle, because neighboring crowns connect together earlier as the zenith angle increases.
Thus, the forest is more similar to a continuous canopy at larger zenith angles.This continuous effect is more apparent for larger crown depth and radius ratios, and it occurs earlier in the solar zenith dimension.

Crown Shape
The tree crown was simplified as an ellipsoid in the proposed model.One convenience of this approach is that different species can be discriminated through the crown shape factor / 2 h R combined with horizontal radius, R.
To highlight crown shape influence on canopy Ω(θ), we kept the crown column constant while changing the crown shape factor (Figure 18).Table 2 shows that crown depth, radius, and their ratio change synchronously.
As the ratio between crown depth and radius increases, there is more obvious increasing Ω(θ) trend with zenith angle.When 2 / 2 h R = , i.e., the crown is spherical, there is only slight Ω(θ) increase with increasing zenith angle (~0.2 in this example).However, when 8 / 2 h R = , implying that a deep canopy, Ω(θ) increased with growing zenith angle, and approaches 1.
At near nadir direction, lower crown depth and radius ratios correspond to larger horizontal projection for each crown when the column is fixed.This larger projection corresponds to lower gap fraction and results in larger Ω(θ).However, curves with higher crown depth and radius ratios demonstrate Ω(θ) increase with growing zenith angle, because neighboring crowns connect together earlier as the zenith angle increases.
Thus, the forest is more similar to a continuous canopy at larger zenith angles.This continuous effect is more apparent for larger crown depth and radius ratios, and it occurs earlier in the solar zenith dimension.

Tree Distribution Pattern
Non-random tree distribution causes clumping and more apparent shadowing among crowns, which also affects canopy Ω(θ) at large zenith angles.Two tree distribution functions were compared: Poisson and Neyman distributions, since these represent random and non-random distributions, respectively [46].
Neyman distribution included m 2 = 1, 2, 5, 10, 20, as shown in Figure 16.When m 2 = 1, the gap fraction distribution is similar to, but not the same as, the Poisson distribution.As m 2 increases, the clumping effect between trees and intra-crown gap fraction increase.However, this combined effect on Ω(θ) is not large.When m 2 ≤ 5, Ω(θ) difference between sample lines <0.02.This demonstrates that the slight non-random distributions in natural forests has little influence on canopy level Ω(θ).But, if the non-random distribution is dominant, its influence on Ω(θ) cannot be ignored.For example, the Ω(θ) difference between the sample line with m 2 = 5 and that with m 2 = 20 can reach 0.08 at large zenith angles.

Trunk Height
The proposed algebraic models for forest clumping index estimation parameterized tree trunk influence.The influence of tree trunk on gap fraction would be great for high trunk tree species, such as the black spruce forest that is considered above.
Figure 17 shows that the influence of trunk is more important at higher zenith angles, because the fraction of visible trunk grows with zenith angle from only the cross section at nadir direction.Higher trunk fraction in the field of view leads to lower gap fraction, causing increased Ω(θ) for constant average LAI.The variation magnitude depends on the forest scenario, which is larger than 0.1 in this example and it should not be neglected.
Note that increased trunk size would increase the equivalent average LAI, since the trunk is one of the vegetation elements, but the average LAI was set as constant in the simulation (Figure 17) to realize monofactor sensitivity analysis and highlight trunk height influence on directional Ω(θ).

Crown Shape
The tree crown was simplified as an ellipsoid in the proposed model.One convenience of this approach is that different species can be discriminated through the crown shape factor h/2R combined with horizontal radius, R.
To highlight crown shape influence on canopy Ω(θ), we kept the crown column constant while changing the crown shape factor (Figure 18).Table 2 shows that crown depth, radius, and their ratio change synchronously.
As the ratio between crown depth and radius increases, there is more obvious increasing Ω(θ) trend with zenith angle.When h/2R = 2, i.e., the crown is spherical, there is only slight Ω(θ) increase with increasing zenith angle (~0.2 in this example).However, when h/2R = 8, implying that a deep canopy, Ω(θ) increased with growing zenith angle, and approaches 1.
At near nadir direction, lower crown depth and radius ratios correspond to larger horizontal projection for each crown when the column is fixed.This larger projection corresponds to lower gap fraction and results in larger Ω(θ).However, curves with higher crown depth and radius ratios demonstrate Ω(θ) increase with growing zenith angle, because neighboring crowns connect together earlier as the zenith angle increases.
Thus, the forest is more similar to a continuous canopy at larger zenith angles.This continuous effect is more apparent for larger crown depth and radius ratios, and it occurs earlier in the solar zenith dimension.

Conclusions
This study investigated the directional clumping effects from typical discontinuous vegetation canopies, including crop and forest.The results verified that the angular variation of the clumping index not only exists, but also reaches notable magnitude for practical cases.
We proposed analytical models to predict the Ω(θ) of crop and forest using structural parameters.The models were validated using simulations for near real cases, and also compared with several sets of ground measurements collected from field experiments and literature.The comparisons provide confidence that the calculated and simulated gap fractions fit ground measurements well.
The analytical model derived from the definition shows that crop Ω(θ) is mainly determined by the row structures, LAI and leaf angle distribution.The effect of these parameters on Ω(θ) will decrease at large view zenith angles.In contrast, forest Ω(θ) is dominated by tree density, tree distribution pattern, crown shape, and trunk size.Their influence on the variation magnitude of Ω(θ) are all not negligible.
Is it necessary to introduce Ω(θ)?In most practical applications of vegetation remote sensing, average clumping index Ω has played the leading role and provided the bridge from effective LAI to actual LAI.Since Ω has integrated the overall hemispherical influence, it could be conveniently rendered as a unified global map.However, the large Ω(θ) variation in the angular dimension and the influence of clumping index on LAI and GPP prompt us to study Ω(θ) systematically.We intend to investigate the proposed approach to further assess the difference between LAI and GPP products that are derived using Ω or Ω(θ), and validate the results using abundant ground measurements.This expansion will be useful to judge whether Ω(θ) should be introduced for standard LAI and GPP production.Another potential option is to apply Ω(θ), as the reference to select an optimal average Ω.
The current research still suffers from many problems and limitations.
(1) The crop model was developed for row structure with homogenous plant elements within the hedgerow.It is not suited for other structures or intercropping of multiple crops.The forest model has not considered the complementary tree species.The tree distribution pattern, LAI, and gap fraction will be quite different in a scenario with the mixture of complementary trees [53].(2) The influence of the crown shape model needs further investigation.Currently, we have assumed the tree crown as vertical ellipsoid.For most species, the crowns can quite well be represented by ellipsoids; however, some exceptions, such as the boreal needle forest, can be fit better by cones McPherson [54] stated that the mean difference between the crown volume measures from the assumptions of crown shape as paraboloids, vertical ellipsoids, and horizontal ellipsoids are 10%.In tropical forest, the multi-stem trees are very common.The current crown shape model might be not valid for them.The further effect of crown shape model on the clumping index needs to be discussed for different species.(3) Further comprehensive sensitivity analysis is expected.The analysis in this paper is restricted to be a mono parameter analysis and the variation range of all the parameters is independent.This is only the first step to observe the influence of each single parameter under specific conditions.In the natural environment, there is a size-density allometry of plants under self-thinning as the resources for growing are limited.The crown projection area scales to stem diameter [53].
Besides, LAI and LAD many vary simultaneously.Therefore, a multiple simultaneous analysis is needed to avoid the risk to misconstrue simple size effects as changes in the crown morphology.
Further research is required to refine this analysis.The analytical directional clumping index models would be applied in the unified model of bidirectional reflectance distribution function for the vegetation canopy [40].With the Unified BRDF model, the leaf area index can be analytically estimated.Then, further validation using ground measurements and inter-comparison with other products will demonstrate whether or how much the introduction of directional clumping index would improve the LAI accuracy.
c and d corresponds to the local peak and valley value.

Figure 1 .
Figure 1.Illustration of effective path length variation with view zenith angles.Figure 1. Illustration of effective path length variation with view zenith angles.

Figure 1 .
Figure 1.Illustration of effective path length variation with view zenith angles.Figure 1. Illustration of effective path length variation with view zenith angles.
c and d corresponds to the local peak and valley value.

Figure 1 .
Figure 1.Illustration of effective path length variation with view zenith angles.

Figure 3 .
Figure 3. Calculation processes for forest average gap fraction.
Set key canopy structure parameters Simulate the tree distribution Assume incident radiation at θ Split the canopy into unit cylinders along θ Calculate the gap fraction within each cylinder Tabulate statistics of the canopy gap fraction

Figure 3 .
Figure 3. Calculation processes for forest average gap fraction.

Figure 6 .
Figure 6.Relative discrepancy between simulated and modeled crop clumping index for different zenith angles at specific azimuth angles.

22 Figure 6 .
Figure 6.Relative discrepancy between simulated and modeled crop clumping index for different zenith angles at specific azimuth angles.

Figure 7 .
Figure 7. Measured and simulated gap fraction from an old black spruce forest, near Candle Lake, Saskatchewan.Parameters are referred to Column 1 in Table2.

Figure 7 .
Figure 7. Measured and simulated gap fraction from an old black spruce forest, near Candle Lake, Saskatchewan.Parameters are referred to Column 1 in Table2.

22 Figure 6 .
Figure 6.Relative discrepancy between simulated and modeled crop clumping index for different zenith angles at specific azimuth angles.

Figure 7 .
Figure 7. Measured and simulated gap fraction from an old black spruce forest, near Candle Lake, Saskatchewan.Parameters are referred to Column 1 in Table2.

Figure 8 .
Figure 8. Simulated distribution pattern of black spruce (Neyman distribution, m 2 = 4, 1000 trees distributed in 4 ha).The white circles represent the nadir projection of the crowns to mark the tree positions, while the black represents the soil background.

Figure 9 .
Figure 9. Simulated and measured gap fraction for Larch stands.Parameters are listed in Column 2 in Table2.

Figure 9 .
Figure 9. Simulated and measured gap fraction for Larch stands.Parameters are listed in Column 2 in Table2.

Figure 15 .
Figure 15.Clumping index variation with zenith angle for different tree density (crown coverage).N = total tree count within 4 ha.Values in brackets represent once-counted crown coverage, which varies according to the actual scenario.Other simulation parameters were as listed in Table 2, column 3.

Figure 15 .
Figure 15.Clumping index variation with zenith angle for different tree density (crown coverage).N = total tree count within 4 ha.Values in brackets represent once-counted crown coverage, which varies according to the actual scenario.Other simulation parameters were as listed in Table 2, column 3.

Figure 16 .
Figure 16.Clumping index variation with zenith angle for different tree distribution patterns.Simulation parameters as shown inTable 2, column 4.

Figure 16 .
Figure 16.Clumping index variation with zenith angle for different tree distribution patterns.Simulation parameters as shown inTable 2, column 4.

Figure 17 .
Figure 17.Clumping index variation with zenith angle for different tree trunk heights.Simulation parameters as shown in Table 2, column 5.

Figure 18 .
Figure 18.Clumping index variation with zenith angle for different crown shape factors.Simulation parameters, as shown inTable 2, column 6.

Figure 18 .
Figure 18.Clumping index variation with zenith angle for different crown shape factors.Simulation parameters, as shown inTable 2, column 6.

Table 1 .
List of symbols used in the clumping index (Ω) models.
θ Mean projection of unit foliage area along the direction of zenith angle θ dimensionless A/A 1 /A 2 /H Row period width, row width, gap width, and crop height, respectively.A = A 1 + A 1 for any ϕ m

Table 2 .
Forest structure parameters used in the current paper.

Table 2 .
Forest structure parameters used in the current paper.
Figure 17.Clumping index variation with zenith angle for different tree trunk heights.Simulation parameters as shown in Table 2, column 5.

Table 2 ,
column 5. Clumping index variation with zenith angle for different tree trunk heights.Simulation parameters as shown in Table2, column 5.