One-dimensional models of radiation transfer in heterogeneous canopies: A review, re-evaluation, and improved model

Despite recent advances in the development of detailed plant radiative transfer models, the large-scale canopy models generally still rely on simplified one-dimensional (1D) radiation models based on assumptions of horizontal homogeneity, including dynamic ecosystem models, crop models, and global circulation models. In an attempt to incorporate the effects of vegetation heterogeneity or “clumping" within these simple models, an empirical clumping factor, commonly denoted by the symbol Ω, is 5 often used to effectively reduce the overall leaf area density/index value that is fed into the model. While the simplicity of this approach makes it attractive, Ω cannot in general be readily estimated for a particular canopy architecture, and instead requires radiation interception data in order to invert for Ω. Numerous simplified geometric models have been previously proposed, but their inherent assumptions are difficult to evaluate due to the challenge of validating heterogeneous canopy models based on field data because of the high uncertainty in radiative flux measurements and geometric inputs. This work provides a critical 10 review of the origin and theory of models for radiation interception in heterogeneous canopies, and an objective comparison of their performance. Rather than evaluating their performance using field data, where uncertainty in the measurement model inputs and outputs can be comparable to the uncertainty in the model itself, the models were evaluated by comparing against simulated data generated by a three-dimensional leaf-resolving model in which the exact inputs are known. A new model is proposed that generalizes existing theory, is shown to perform very well across a wide range of canopy types and ground cover 15

Ω in the exponential argument of Beer's law that effectively scales the radiation attenuation coefficient based on the level of vegetation clumping (Nilson, 1971;Chen and Black, 1991;Black et al., 1991). Clumping usually results in an over estimation of radiation absorption when a model based on Beer's law is used (Ponce de León and Bailey, 2019), and thus setting Ω < 1 reduces the effective attenuation coefficient within Beer's law, which corrects for this over estimation.
Despite the wealth of available models for quantifying radiative transfer in heterogeneous canopies, a critical knowledge gap 55 still exists in which it is usually unclear which model is suited for a particular application, and even a general sense of the errors associated with certain model assumptions is often unknown. Models are commonly selected for historical reasons, based on ease of implementation, availability of computational resources versus domain size, or on perceived errors given the particular model assumptions. This uncertainty is driven by the fact that obtaining robust validation data is exceptionally difficult, and often uncertainty in model inputs is comparable to uncertainty in the model itself. Offsetting errors and coefficient "tuning" can 60 lead to models that perform exceptionally well in a particular case, but may produce unacceptably large errors when applied generally. These difficulties have led to a number of model inter-comparison exercises in which simulations are performed of synthetic or artificial canopy cases (Pinty et al., 2001(Pinty et al., , 2004Widlowski et al., 2007Widlowski et al., , 2013. This eliminates ambiguity in model inputs to enable a more objective comparison, however, the "exact" solution is still unknown. This paper presents a critical re-evaluation of the theoretical basis of simplified one-dimensional (1D) models of radiation transfer in heterogeneous canopy based on Beer's law. Due to the difficulties in objectively comparing and evaluating models based on field data, the performance of various models was explored by applying them in virtually-generated canopies where the inputs are exactly known, and comparing against the output of a detailed leaf-resolving model. The goal of the study was to better understand the implications of radiation model assumptions and uncertainty in model inputs in a wide range of canopy geometries in order to guide model selection in future applications.  75 where I(r; s ) is the radiative intensity at position r along the direction of propagation s , κ and σ s are the radiation absorption and scattering coefficients of the medium, respectively, I b (r) is the blackbody intensity of emission at position r, Φ(r; s , s) is the scattering phase function at position r for incident direction s and scattered direction s, and dΩ is a differential solid angle.
This equation can be integrated along s from a distance of r = 0 to r to yield I(r; s ) = I(0; s )exp (−κr) , The attenuation coefficient can be interpreted physically as the cross-sectional area of radiation-absorbing objects projected in the direction s per unit volume of the medium. For a medium of leaves, the attenuation coefficient is given by where a is the one-sided leaf area density (m 2 leaf area per m 3 canopy), and G(s ) is the fraction of total leaf area projected in the direction of s . It is frequently assumed that the attenuation coefficient does not have azimuthal dependence, and therefore the G-function can be written as G(θ), where θ is the zenithal angle of the direction of radiation propagation. The probability of a beam of radiation, inclined at an angle of θ, intersecting a leaf within a homogeneous volume of vegetation after propagating 90 a distance of r can thus be written as For a canopy that extends indefinitely in the horizontal, the propagation distance r can be re-written in terms of the canopy height h as r = h/cos θ. Substituting this relation for r, and noting that the leaf area index (LAI) is defined as a h = L (if a is constant), yields the common form of Beer's law applied to plant canopy systems

Application of Beer's law in heterogeneous canopies
As introduced previously, Eq. 6 is only valid within a homogeneous volume of vegetation, i.e., leaves are uniformly distributed in space. In reality, essentially all canopies have heterogeneity or "clumping" at a number of scales. There is inherent clumping at the leaf scale, as leaves are discrete surfaces unlike arbitrarily small gas molecules, and thus even at this scale the assumptions 100 of Beer's law are violated. Above the leaf scale, leaves are grouped around shoots or spurs, creating heterogeneity at a larger scale. Shoots are grouped around main branches or the plant stem, creating yet another scale of heterogeneity. Large-scale clumping typically exists due to gaps between individual trees, or clearings in the canopy.
In a strict sense, each of these scales of heterogeneity violates the assumptions of Beer's law. One obvious means of dealing with this heterogeneity is to use a more complicated model that explicitly resolves the important scales of heterogeneity, such 105 as a "multilayer" model that resolves heterogeneity in the vertical direction (e.g., Meyers and Paw U, 1987;Leuning et al., 1995), a 3D model that resolves plant-scale heterogeneity (e.g., Wang and Jarvis, 1990;Cescatti, 1997;Stadt and Lieffers, 2000), a 3D voxel-based model that resolves sub-plant heterogeneity (e.g., Kimes and Kirchner, 1982;Sinoquet et al., 2001;Gastellu-Etchegorry et al., 2004;Bailey et al., 2014), or a 3D leaf-level model that resolves heterogeneity at the leaf scale (e.g., Pearcy and Yang, 1996;Chelle and Andrieu, 1998;Bailey, 2018;Henke and Buck-Sorlin, 2018). However, each of these 110 models incurs a level of computational expense that may be unacceptable for global-scale models or crop models, which require high efficiency that is provided by simpler models based on Beer's law. As a compromise, numerous geometric models have been proposed that calculate radiation interception using analytical geometric solutions along with simplifying assumptions of the basic shape of canopy elements, such as a hedgerow crop (Tsubo and Walker, 2002;Annandale et al., 2004), spherical/ellipsoidal crowns (Norman and Welles, 1983;Li and Strahler, 1988;Ni-Meister et al., 2010), or conical crowns 115 (Kuuluvainen and Pukkala, 1987;Van Gerwen et al., 1987;Li and Strahler, 1988). However, these models do not necessarily generalize to an arbitrary canopy, and in many cases the choice of the average geometric input parameters may be unclear and thus their determination may require an empirical inversion (e.g., Li and Strahler, 1988).

"Ω" canopy clumping factor approach for incorporating vegetative heterogeneity
The issue of incorporating the effects of clumping in Beer's law models gained a heightened level of attention in the early 120 1990s from investigators looking to use radiation measurements to invert Beer's law for leaf area index (LAI) values Black et al., 1991;Chen and Black, 1993). Canopy non-randomness or clumping causes underestimation of LAI values inferred in this way, and thus there was a pressing need for a theoretical formulation that could remove the effects of clumping within the inversion procedure, which was also mathematically simple enough that it could be easily inverted for LAI. In an early attempt at applying such a correction, Chen and Black (1991) introduced an empirical coefficient within Beer's 125 law which was termed the "clumping index" and denoted by Ω: Chen and Black (1991) references the early work of Nilson (1971), who originally derived this relationship (his Eq. 25, with the "clumping" parameter denoted by λ). Nilson (1971) derived this relationship for "stands with a clumped dispersion of foliage" using a Markov chain model. The assumption was that the probability of interception at any point also depends to some degree 130 on the probability of interception at a previous location, with the level of dependence described by Ω (or λ using the notation of Nilson (1971)). This approach, however, does not explicitly allow for large-scale gaps in vegetation such as crown-scale clumping, but rather assumes that there is some quasi-homogeneous medium with regular and repeating variation in vegetation density.
Despite the seemingly attractive simplicity of the clumping factor approach, it has many limitations. Foremost of these limitations is that in general the clumping factor Ω has a strong dependence on every other variable in Eq. 7, namely θ, L, and G (Chen et al., 2008;Ni-Meister et al., 2010). In effect, specifying Ω requires knowledge of P . Thus, if P is already known, there is no need to determine Ω and apply Eq. 7 in the first place. Empirical modeling of Ω (e.g., Kucharik et al., 1997;145 Campbell and Norman, 1998;Kucharik et al., 1999) is effectively a matter of empirically modeling P . If the leaf orientation and heterogeneity are both isotropic, then G and Ω would no longer have θ dependence. In this case, the product G Ω L becomes inseparable and the canopy begins to appear homogeneous with attenuation determined by the value of G Ω L. The probability of interception can then be written as

150
where Ω = G Ω L = const. Thus, if P is known for any particular θ, Ω can be determined.

Direct (collimated) radiation component
Beer's law is only explicitly valid in a medium in which the probability of radiation interception is homogeneous over some discrete scale. Thus, in order to apply Beer's law in heterogeneous vegetation, we must segment the canopy into sections 155 over which we can assume that the vegetative elements are homogeneous in space. In a typical canopy, there may be multiple scales over which this is applicable ( Fig. 1). At the leaf scale, the probability of interception is approximately homogeneous, and is equal to the leaf absorptivity. For a "clump" of leaves (e.g., a shoot), the probability of interception may be assumed approximately homogeneous, and is given by Beer's law (Eq. 6). There may be large gaps between individual shoots, which invalidate Eq. 6, but the overall distribution of shoots across a crown may be approximately homogeneous, and thus the 160 probability that a photon intersects an individual shoot may follow Beer's law but with an augmented attenuation coefficient. Finally, there may be large gaps between crowns that invalidate Eq. 6, but this heterogeneity may be regular and repeating, and thus the probability that a photon intersects an individual crown may be described by Beer's law.
The cumulative probability of photon interception can be viewed as an aggregation of many "clumps" consisting of a homogeneous medium of elements at a smaller scale. For each clump, the probability of a photon intersection is assumed homo-165 geneous in space, and thus the probability of interception can be assumed constant. In this case, the cumulative probability of interception over all clumping levels is where P is the cumulative probability of interception over all scales, N c is the total number of clumping levels, and P i is the probability of intersecting the i th clumping level.

170
The product of Eq. 9 could be applied in any number of ways depending on the scale at which various probabilities of interception are known. In the example given below, we will follow an approach similar to Nilson (1999) in which the probability of interception within a single plant crown is determined, then repeated N c times for each crown in the canopy that the beam of radiation traverses. Accordingly, the canopy is segmented into crown "envelopes" (cf. Nilson, 1992), each of which is conceptualized to a volume encompassing all vegetation within the plant, inside which it is assumed that vegetation 175 is homogeneous.
The probability of a beam of radiation intersecting a single crown is the product of the probability of intersecting the crown envelope and the probability of intersecting a leaf within the envelope. At a solar zenith angle of zero, the probability of intersecting the crown envelope is given by the ground cover fraction f c , which is the area of the crown envelope shadow at a solar azimuth of zero S(0) divided by the average plan area of ground associated with a single plant. For spherical or 180 cylindrical crowns of radius R and spacing s, the ground cover fraction is f c = S(0)/s 2 = πR 2 /s 2 . The probability that a beam intersects a leaf within a single crown is simply given by Eq. 5, with r being the path length through the crown. Because of the non-linearity of this equation, we cannot simply compute the average r over the crown and substitute it into Eq. 5.
Rather, Eq. 5 should be weighted by the probability that a beam path length through the crown is equal to some value r. Thus, the probability P that a ray passing through a crown envelope intersects a leaf can be written as where p(r; θ) is the probability that the beam path length through a crown is equal to r for a solar zenith angle of θ. Li and Strahler (1988) showed that for a sphere of radius R, p(r) = r/2R 2 (no θ dependence for spheres), with 0 ≤ r ≤ 2R and thus the limits on the integral in Eq. 10 become 0 to 2R. For other shapes, analytical expressions for p(r; θ) become tedious or impossible to derive. However, the integral expressions given by Li and Strahler (1988) for p(r; θ) can be evaluated numerically.

190
The approach used here was to perform line-cylinder intersection tests (which are analytical Suffern, 2007), for a large number of lines inclined in the direction of the sun, which allows for population of a probability distribution. A similar approach could be used for ellipsoids (equations given in Norman and Welles, 1983).
In order to use the above expression to calculate the probability of interception for an entire canopy of repeated crowns, we must assume a statistical distribution that describes the probability of intersecting a crown envelope within a canopy. The 195 two distributions that are commonly chosen are the binomial distribution (Nilson, 1971(Nilson, , 1999 or Poisson distribution (Nilson, 1971;Li and Strahler, 1988;Nilson, 1999). If a (positive) binomial distribution is assumed, the interception probability for the entire canopy is where N c is the number of crowns intersected by the radiation beam. Effectively this amounts to saying that the probability of

200
NOT intercepting an individual crown is 1 − f c P , which is compounded N c times. When the solar zenith is zero, this means that N c = 1 and thus Eq. 11 correctly yields a probability of intersection of f c P .
If a Poisson distribution is assumed, the interception probability for the entire canopy is Application of this model effectively assumes that crowns act like a homogeneous medium of objects that repeats in all direc-205 tions, analogous to individual molecules in a gas. When the sun direction is near zenith, this assumption is poor as the canopy consists of only a single layer of objects. In this case, the probability of intersection should be f c P , but Eq. 12 incorrectly gives a value of 1 − exp (−f c P ), which is always less than f c P .
This approach for estimating N c works well if the crowns are randomly or uniformly distributed in space, and thus N c is azimuthally symmetric. If crowns are oriented in rows, the relationship for N c changes with azimuth. To account for this, we propose the following. Consider the case in which crowns are spaced at a distance of s p (plant spacing) in the direction of 215 radiation propagation, and spaced at a distance of s r (row spacing) in the direction normal to the direction of propagation.
The azimuthally symmetric model of N c = S(θ)/S(0) can be applied, provided that the asymmetry in the actual ground cover fraction is properly accounted for. This is accomplished by 1) calculating the ground cover fraction using the crown spacing in the direction of radiation propagation, which in this example is f c = S(0)/s 2 p , and 2) multiplying the final interception probability by the ratio of the isotropic crown footprint area (s 2 p in this example) to the actual crown footprint area (s p s r in this 220 example).
To generalize this approach, we take the azimuthally symmetric plant spacing s to be s = s r sin 2 ϕ + s p cos 2 ϕ, where ϕ is the azimuthal angle between the sun direction and the row direction. Equation 11 can then be generalized to It is noted that when s r = s p , Eq. 13 reduces back to Eq. 11. Table 1 provides a summary of inputs and equations needed to 225 implement Eq. 13 for different geometries.

Diffuse radiation component
As introduced above, the RTE (Eq. 1) and thus Beer's law (Eq. 6) are only explicitly valid along a single direction of radiation propagation, and therefore these equations as written can only be applied for collimated radiation (e.g., direct solar radiation).
It is common to adapt Beer's law for diffuse radiation conditions by substituting a modified diffuse radiation attenuation 230 coefficient that is usually assumed to be constant for a particular canopy (e.g., DePury and Farquhar, 1997; Wang and Leuning, 1998;Drewry et al., 2010). However, for the reasons above, this approach is not consistent with the assumptions of Beer's law.
A more robust approach for representing incoming diffuse radiation from the sky is to apply Beer's law to any given direction in the sky and integrate across the upper hemisphere according to 235  (0) s Adjusted effective spacing for canopies with rows oriented at an angle of ϕ relative to the sun s r sin 2 ϕ + s p cos 2 ϕ** P Probability of intersecting a leaf within a crown **If plant spacing is azimuthally symmetric (e.g., random), s = sr = sp. Beer's law with empirical clumping factor model (see Campbell and Norman, 1998) HOM Eq. 6 Standard Beer's law -leaves randomly positioned in space where P diff is the fraction of incoming diffuse radiation intercepted by the canopy, P (θ, φ) is the fraction of radiation intercepted by the canopy for radiation originating from the spherical direction of (θ, φ) (such as that calculated by Eq. 13), and f d (θ, φ) is a weighting factor to account for anisotropic incoming diffuse radiation. If the canopy is assumed azimuthally symmetric, this equation reduces to Thus, calculation of diffuse radiation interception is simply a weighted average of collimated radiation originating from the upper hemisphere. In the case of a uniform overcast sky, f d = 1 and thus the amount of collimated radiation originating from some direction θ is weighted by cos θ sin θ.

Specification of model inputs
For essentially all of the models considered in this work, the model input parameters are 1) the leaf G-function, 2) leaf area index or density, 3) the relative density of crowns, and 4) a mathematical description of the crown envelope. Most commonly, the crown envelope is assumed to be spherical or ellipsoidal, and thus it is described by its radii. Objectively specifying the crown envelope is generally more of a challenge than it may seem, particularly when recalling that we should define the 250 envelope such that it can be assumed that the vegetation inside the envelope is uniformly distributed in space. Typical crowns have shoots and branches that create an additional scale of clumping, and create irregularly shaped crowns. Figure 2 shows an example of potential spherical and ellipsoidal crown envelopes that could be chosen for a few trees. Clearly, no matter how the envelope is defined, there is a significant fraction of the envelope that contains open spaces with no leaves, which violates our model assumptions. More complicated envelopes can be derived that better fit the shape of the crown (Nilson, 1999), but these 255 cases require numerical integration in order to calculate the relevant model inputs, namely S(θ) and p(r).
3 Test case set-up

Overview
While it is possible to test the above modelling framework using field data, this approach is severely limited by the lack of systematic variation in canopy architecture as well as experimental errors that can become convoluted with model errors. As 260 such, it can be difficult if not impossible to use field data to rigorously evaluate and diagnose issues within models. In this work, an alternative approach was used to evaluate model performance, which was based on the use of a detailed 3D leaf-resolving model to simulate radiation absorption in virtually-generated canopies. The advantage of this approach is that arbitrary canopy geometries can be generated in which the exact geometry is known, which provides the necessary inputs for a 1D model.
The limitation of course is that results are confined within the assumptions and accuracy inherent in the chosen 3D model.

265
In order to minimize this limitation, the ray-tracing-based model of Bailey (2018) was used as implemented in the Helios 3D modelling framework (Bailey, 2019). Provided that simulated surfaces are isotropic absorbers, Bailey (2018) showed that this model converges exponentially toward the exact solution as the number of rays is increased. It has also been shown to converge to the solution given by Beer's law for a truly homogeneous canopy (Ponce de León and Bailey, 2019). Since it was verified that further increasing the number of rays did not significantly affect results, we considered the 3D model solution to be the 270 reference or "exact" solution against which the various 1D models could be compared. The canopy-level intercepted radiation flux was calculated from the 3D model output as where R ↓ is the above-canopy solar radiation flux on a horizontal surface, A g is the horizontal area of the canopy footprint, R i is the radiative flux incident on the i th of N p total vegetative elements, and A i is the one-sided surface area of the i th vegetative 275 element.
A number of test cases were formulated to progressively test different aspects of each of the models given in Table 2. For each of the test cases below, all surfaces were black, and the ambient diffuse radiation flux was set to zero in order to maintain consistency with these non-geometric assumptions inherent in Beer's law. The above-canopy solar radiation flux was modelled using the REST-2 model (Gueymard, 2003). The model was run at an hourly timestep for Julian day 79, and the position on the 280 earth was the equator. This date and location were chosen because a single diurnal cycle provides a full range of solar zenith angles ranging from 0 to π/2. This means that it is not explicitly necessary to evaluate performance under diffuse conditions because, as was illustrated by Eq. 15, the diffuse radiation flux is simply a weighted average of the flux originating from all zenithal directions.
Agreement between each of the 1D models and the 3D model were analyzed graphically, and quantitatively using the index 285 of agreement (Willmott, 1981), which is defined mathematically as In order to isolate the effects of crown-scale clumping, a test case was considered in which the canopy consisted of solid, opaque spheres of radius R = 5 m with varying spacing (Fig. 3a). In the first configuration, spheres were arranged randomly with an average spacing in the horizontal direction of s/R = 2, 3, 4, and 6. It should be noted that the placement of spheres was (a) not truly random, as crowns were not allowed to overlap. The second configuration placed the spheres in a non-random row orientation in which the plant spacing in the row-parallel direction was s p /R = 2, 3, 4, and 6 and the row spacing was s r /R = 295 4, 6, 8, and 12. Rows were oriented in either the north-south or east-west directions, and crowns were not allowed to overlap.
The spheres were positioned on the ground surface, and thus the canopy height was h = 2R. The 1D models were evaluated with L = ∞ since the spheres were solid. Only the models BINOM, NIL99_B, and NIL99_P were considered for this case because they separately account for crown and sub-crown intersection.

300
To test generalization to geometries with anisotropic crowns, a case was considered with a canopy of solid, opaque cylinders of radius R = 5 m and height H = 2R (Fig. 3b). The set-up was essentially the same as Case #1, with random, north-south, and east-west arrangements of crowns with the same set of spacings.

Case #3: Canopy of uniformly distributed leaves in spherical crowns
In order to test the combined effects of crown-scale clumping and leaf-scale attenuation, a test case was considered in which the 305 canopy consisted of spherical crowns containing homogeneous and isotropic vegetation elements (Fig. 3c). Spherical crowns of radius R = 5 m were generated in which rectangular leaves of size 0.1×0.1 m 2 were arranged randomly inside the crown with a uniform spatial distribution, and were randomly oriented following a spherical distribution (G = 0.5). For brevity, only the random crown arrangement is presented, which had the same average spacing as in Case #1 above. The leaf area density within each crown was set at a = 0.5 m −1 , and the whole-canopy LAI can be calculated as L = 4πR 3 a 3s 2 , which gives LAI 310 values ranging from 0.3 to 2.6.

Case #4: Canopy of uniformly distributed leaves in cylindrical crowns
Similar to Case #3, an additional case was considered consisting of cylindrical crowns of radius R = 5 m and H = 2R, filled with uniformly distributed leaves (Fig. 3d). All other parameters are the same as Case #3. For cylindrical crowns, the wholecanopy LAI is L = πR 2 Ha s 2 , which gives LAI values ranging from 0.44 to 3.9.

Case #5: Canopy of trees
In order to test the models for more realistic canopy architectures, canopies of trees were constructed using the procedural tree generator of Weber and Penn (1995) (Fig. 3e) as implemented in the Helios 3D modelling framework (Bailey, 2019). The trunk and branches consisted of a triangular mesh of elements forming cylinders, and leaves were rectangles of size 0.12×0.03 m −2 that were masked into the shape of a leaf using the transparency channel of a PNG image of a leaf (see Bailey, 2019). The tree 320 crown envelope was approximately spherical in shape, but the spatial distribution of leaves was non-uniform. Leaf angles were sampled from a uniform distribution, and thus G = 0.5. Note that in calculating radiation interception, a distinction was not made between branches or leaves, but rather total attenuation was used. The tree height was approximately h = 6.5 m, and trees were randomly arranged with average spacing in the horizontal direction of s = 5, 7.5, 10, 15, and 20 m. The canopy leaf area index values were L = 3.86, 1.57, 0.97, 0.39, and 0.24. An effective crown radius was estimated to be R = 2.9 m, which was 325 the value that gave the best predictions of radiation interception at a solar zenith angle of zero. Figure 3e shows a visualization of the assumed crown envelope based on a best fit of both the BINOM and NIL99_BINOM models to the exact intercepted flux. The leaf area density was variable in space, and thus an effective density within the crown was substituted into the 1D equations, which was calculated as a = 3Ls 2 4πR 3 .

330
A potato plant canopy was generated to test the models in more realistic non-tree canopies (Fig. 3f). Similar to the tree canopies, the stem consisted of a mesh of triangles, and leaves were texture-masked rectangles. The crown envelope could be considered roughly cylindrical, with leaves of variable size distributed non-uniformly within the cylinder. The leaf angle distribution was highly anisotropic, with G ranging from 0.87 when the sun direction was vertical to 0.23 when the sun direction was horizontal.
The effective dimensions of the crown envelope were estimated to be R = 0.5 m and H = 0.75 m, which is visualized in Fig. 3f.

345
For the east-west row-oriented configuration (Fig. 4e-h), performance of the BINOM model was nearly the same as in the randomly oriented case. As would be expected, the performance of the NIL99_B model decreased in the east-west row-oriented case as it assumes an azimuthally symmetric distribution of crowns (e.g., d decreased from about 1.0 to 0.94 in the sparsest case). The performance of the Poisson model (NIL99_P) actually increased slightly for the east-west orientation, however this is likely the result of offsetting errors. As evidenced by the NIL99_B results, an east-west row orientation appeared to 350 cause an over prediction of the intercepted flux, which offsets some of the under prediction due to the assumption of a Poisson distribution in the probability of crown intersection.
Overall, the BINOM model performed equal to or better than the NIL99_B and NIL99_P for every canopy configuration.
The lowest d value for the BINOM, NIL99_B, and NIL99_P models, respectively, were 0.98, 0.94, and 0.81. It is noted that for   Figure 4. Flux of intercepted radiation versus solar zenith angle for a half diurnal cycle in a canopy of solid spheres (see Fig. 3a). Varying ground cover fractions fc were simulated (rows in figure) for three different plant arrangements: randomly spaced (a-d), east-west row orientation (e-h), and north-south row orientation (i-l). The "exact" flux is given by the output of the 3D model ( ), which is compared against the 1D models BINOM ( ), NIL99_B (•), and NIL99_P (•).
the case of a canopy of solid objects with random spacing, the BINOM and NIL99_B models are mathematically equivalent, 355 which is confirmed by the results.
4.2 Case #2: Canopy of solid cylinders Figure 5 gives timeseries of the "exact" intercepted radiation flux compared against the simplified 1D models for each of the canopies of solid cylinders with four different densities and three plant arrangements. Trends in model performance were similar as in the solid sphere case (Fig. 4), except that overall performance for all models was decreased slightly. the radiation path length is constant across the crown cross-section, which caused a slight increase in the intercepted flux for NIL99_P (as was the case for the NIL99_B versus BINOM models).
The OM_VAR model had very large errors as the canopy became increasingly sparse. For zenith angles near 90 • , this model had an over prediction as large as the homogeneous model, which decreased toward that of the other models as zenith angle decreased. It should be noted that the OM_VAR is forced to match the exact flux perfectly at θ = 0, and it only needs to model 375 the flux for θ > 0. The OM_CON model, which also is forced to perfectly match the exact flux at θ = 0, performed as well as the BINOM model (d ≥ 0.99). This model, which assumes a constant impact of heterogeneity for θ > 0, performed much better than the OM_VAR model. Such a result is to be expected, given that for spherical crowns, the crown envelope as well as leaf orientation are isotropic. Thus, since Ω(θ) is dependent on both the impact of heterogeneity and G(θ), the fact that both G and the heterogeneity are approximately constant means that Ω is approximately constant.  Figure 5. Flux of intercepted radiation versus solar zenith angle for a half diurnal cycle in a canopy of solid cylinders (see Fig. 3b). Varying ground cover fractions fc were simulated (rows in figure) for three different plant arrangements: randomly spaced (a-d), east-west row orientation (e-h), and north-south row orientation (i-l). The "exact" flux is given by the output of the 3D model ( ), which is compared against the 1D models BINOM the sparsest case, the Poisson model NI10_P still significantly under predicted the flux, whereas the NIL99_P over predicted the flux.
The OM_VAR model had a very large over prediction of radiation interception at high zenith angles, as in the spherical 390 crown case. Unlike in the spherical crown case, the OM_CON model did not perform well, particularly as the canopy became increasingly sparse. When crowns are cylindrical, heterogeneity is no longer isotropic especially as the plant spacing becomes large, and as such interception varies irregularly with θ. As a result, Ω has strong θ dependence and cannot be assumed constant. Figure 8 gives timeseries of the "exact" intercepted radiation flux compared against the simplified 1D models for each of the 395 canopies of randomly spaced trees with four different average densities. Results were quite similar to that of Case #3 (spherical crowns). Model agreement with the reference intercepted flux was reduced slightly overall for each of the models, which is likely due to the fact that the assumptions of spherical crown shape and within-crown homogeneity were not exactly satisfied.

Case #5: Canopy of trees
However, this reduction in performance was small, and the BINOM model still performed very well. Thus, this test confirms generalizability to more realistic tree architectures in which vegetation within the crown envelope is only approximately 400 homogeneous and the crown shape is only approximately spherical.
4.6 Case #6: Non-tree canopy canopy with cylindrical crowns (Case #4), with only slightly reduced overall performance in comparison to Case #4. It is noted 405 that this was the only case with an anisotropic leaf angle distribution, but the high anisotropy in G did not seem to significantly affect model performance.

Discussion
The simplest models, based on an Ω clumping factor that effectively scales the attenuation coefficient, had mixed success depending on the particular case. Assuming a constant Ω factor (OM_CON model) worked quite well in the case that the 410 heterogeneity (crowns) was roughly isotropic. As mentioned previously, Ω encapsulates the effects of G(θ), heterogeneity, and the apparent LAI that results from the heterogeneity to form an inseparable product G Ω L. If this product is isotropic, then it is reasonable to assume a constant Ω. The problem, however, is still that this constant Ω is not known a priori and must be determined from radiation measurements at a particular solar zenith angle (preferably near θ = 0) in order to invert for the appropriate attenuation coefficient.

415
The OM_VAR model (Campbell and Norman, 1998) attempts to model the effects of anisotropic clumping for θ > 0. When the canopy was fairly dense, the OM_VAR model worked fairly well. However, the OM_CON and HOM models generally performed as well or better than the OM_VAR model for those cases, and thus there is no reason to use a variable Ω factor.
As the canopy became increasingly sparse, the performance of the OM_VAR model declined significantly and over predicted interception as the solar zenith angle increased (note that this model also forces the predicted flux to match the reference flux 420 exactly at θ = 0). Kucharik et al. (1999) evaluated the framework behind the OM_VAR model in a number of different canopies and found that it was able to fit the data well, but that the model coefficients were highly species-specific. One issue with their validation approach, which is symptomatic of many field validation studies of heterogeneous canopy radiation models, is that the data was collected in relatively dense canopies where heterogeneity is fairly low overall. The canopies studied in Kucharik et al. (1999) had a ground cover fraction approximately in the range of f c = 0.5 − 0.75, which would place them on the denser 425 end of the canopy cases considered in the present study, which is where the OM_VAR model worked well. However, it was shown that for these cases a constant Ω factor model (OM_CON) or even the homogeneous model (HOM) also worked well.
The assumption that the probability of intersecting a crown envelope follows a Poisson distribution did not work well unless the canopy was very sparse. For most canopy densities, the Poisson models significantly under predicted the absorbed flux. It appeared that when the mean free path of radiation propagation (which is related to the crown spacing) was not significantly (2010) were quite dense with high LAI and ground cover fraction. As illustrated previously, specification of the crown radius can be ambiguous for real trees that are irregularly shaped. If the crown radius is specified based on an envelope encapsulating all branches, this effectively inflates R and f c , which would presumably result in an over prediction of the intercepted flux.

440
However, this over prediction could be offset by applying a Poisson model, which was shown to cause under prediction. It is possible that R was inflated in Yang et al. (2010) (which considered only dense canopy cases), and offsetting errors associated with the Poisson assumption resulted in artificially improved model performance due to offsetting errors. The crown radius in Case #2 is exactly known, so specification of R in that case is not a potential source of model error. In another study, the NI10_P model was compared against other models for a set of virtual canopies where the exact geometry was known (but the 445 exact radiation interception was not), which suggested that the NI10_P model tended to predict higher canopy transmission (lower interception) than the other models (Widlowski et al., 2013), which is also consistent with the results of the present study.
The assumption that crown intersection followed a binomial distribution appeared to hold for all canopy cases considered in this work. The binomial model predicts the correct interception at θ = 0, which corresponds to a single crown layer (N c = 1) 450 and P = 1 − f c , rather than P = 1 − exp(−f c ) in the case of the Poisson model. The BINOM model outperformed all other models for every test considered. The primary difference between the formulation of the BINOM and NIL99_B is that 1) the NIL99_B model assumes that the path length of radiation through an individual crown is constant whereas the BINOM model accounts for variable path length, and 2) the NIL99_B model assumes that crowns are randomly positioned in space and thus crown intersection is azimuthally symmetric, whereas the BINOM model accounts for asymmetry. The assumption of constant 455 path lengths created modest errors, as evidenced by the results of Cases #3 and #4. The assumption of random positioning of crowns in a row-oriented canopy had the potential to create very large errors, as evidenced by the results of Cases #1 and #2.
These errors are caused by the fact that the effective total path length through vegetation in a row-oriented canopy can change significantly with azimuth.
The scope of the results of this study are clearly limited to cases of no scattering, and no diffuse radiation. These impacts 460 were excluded from the study to focus on cases where, aside from heterogeneity in the geometry, the assumptions of Beer's law should be exactly satisfied. Although Beer's law is only valid along a single direction of radiation propagation, and its derivation requires the removal of scattering terms in the RTE, variations have been derived that approximate the effects of scattering and diffuse radiation within a 1D model (e.g., Lemeur and Blad, 1974;Goudriaan and Van Laar, 1994).
It appears likely that many crop models, global ecosystem models, and land surface models over estimate radiation intercep-465 tion by applying the homogeneous Beer's law in heterogeneous environments, which is sure to have important consequences for large-scale flux estimates. Incorporation of the results of this work within these models is straightforward, and requires specification of either the ground cover fraction f c or the planting density and effective crown envelope. Algorithms are readily available for separation of the ground surface and vegetation within aerial images in order to calculate the ground cover fraction (e.g., Gougeon, 1995;Luscier et al., 2006;Laliberte et al., 2007). The results of Cases #5 and #6 suggested that rough 470 estimations of the crown envelope dimensions based on visual inspection could yield reasonable results.

Conclusion
Simplified models of radiation interception in heterogeneous canopies can be readily derived by separating the canopy into hierarchical scales of "clumping" over which the probability of interception can be assumed homogeneous in space over some discrete volume. The results of this work demonstrated that very good predictions of whole-canopy interception can be 475 achieved using simple geometric models that consider only crown-scale and leaf-scale clumping (in the absence of scattering).
The probability of intersecting a plant crown was well represented by a (positive) binomial distribution. This model calculates the probability of not intersecting a leaf within a single crown, and compounds this probability N c times, where N c is the number of crowns a given beam of radiation traverses on its path from the top to the bottom of the canopy. The Poisson models for crown intersection did not perform well unless the canopy was fairly dense, but in this case the effects of heterogeneity are 480 less important and the homogeneous Beer's law also performs well. The results of the model evaluation exercise confirms that the binomial model given in Eq. 13 (BINOM) is the preferable model in all cases considered herein. Inputs to the model can be specified based on measurements of plant geometry, or through inversion if radiation interception measurements are available for a particular solar zenith angle.
Code and data availability. Helios code version 1.0.14 along with associated project files and output files can be downloaded For completeness, model equations taken from the literature are provided below as they were implemented and using the notation adopted in this paper.

A1 Model of Nilson (1999): NIL99_B and NIL99_P
Assuming that the probability of intersecting a crown envelope follows a (positive) binomial distribution, Nilson (1999) gives the probability of intersection for a canopy to be 500 P = 1 − 1 − (1 − P 1 )S(0)/s 2 Nc , where N = S(θ)/S(0), S(θ) is the area of the crown envelope shadow, S(0) is the area of the crown envelope shadow at a solar zenith angle of zero, s is the mean plant spacing, and P 1 is the probability that a beam of radiation does not intersect a leaf within a single crown and is given by If the probability of intersecting a crown envelope is instead assumed to follow a Poisson distribution, Nilson (1999) gives the probability of intersection for a canopy to be P = 1 − exp − S(0) s 2 (1 − P 1 )N c = 1 − exp − S(θ) s 2 (1 − P 1 ) , There are two notable differences between Eq. 13 and the model of Nilson (1999). The first was mentioned above, which is that the expression for N c in Nilson (1999) assumes a random or uniform spatial distribution of crowns, whereas Eq. 13 has 510 been generalized to include row-oriented crowns. Second is that Nilson (1999) assumed a spatially constant within-crown path length in calculating P 1 , whereas P accounts for variable beam path lengths through crowns (Eq. 10).
A2 Model of Campbell and Norman (1998): OM_VAR Kucharik et al. (1997) originally suggested a simple empirical model describing the θ-dependence of Ω, provided that the value of Ω at θ = 0 is known: 515 Ω(θ) = Ω(0) Ω(0) + (1 − Ω(0)exp (−kθ p )) , where Ω(0) is the value of Ω when the solar zenith angle is zero, and k and p are geometric coefficients with p given by (Campbell and Norman, 1998) p = 3.8 − 0.46D, 1 ≤ p ≤ 3.34 and D is the ratio of the crown depth to the crown diameter. The coefficient k is taken to be equal to 2.2 as suggested by 520 Campbell and Norman (1998).
An obvious limitation of this approach is that the value of Ω(0) must be estimated, which usually requires midday solar interception data. It does, however, make modeling easier since this approach guarantees that the correct radiation interception will be predicted around midday.  Li and Strahler (1988). Their approach used geometry of spheres to back-calculate the appropriate clumping factor Ω, which is given by Ω = 3 4τ 0 r 1 − (1 − (2τ 0 r + 1)exp (−2τ 0 r)) 2(τ 0 r) 2 , where 530 τ 0 r = 3G(θ)L 4λπR 2 , and recalling that G is the fraction of leaf area projected in the direction of radiation propagation, L is the canopy leaf area index, λ is the number of plants per unit ground area, and R is the crown radius.
Since Ni-Meister et al. (2010) formulated their model in terms of an Ω clumping factor, the end equations for the models of Ni-Meister et al. (2010) and Nilson (1999) look quite different. However, mathematically they are nearly the same. The primary 535 difference in the formulation is that Nilson (1999)