# Light interception

A conceptual model of light interception which meets the objectives of computational efficiency required by the general model has been developed. The model of light interception is based on the ray-tracing approach (Brunner 1998;Deutschman, Levin, & Pacala 1999), but incorporates Monte Carlo integration and importance sampling to reduce the number of calculations necessary for each call to the light interception model. The model is based on the assumption that the light environment experienced by a point of coordinates [x,y,z] within a forest stand is a function of the light intercepted by the crowns along the rays of azimuthal direction a and zenithal direction z originating from the point. In absence of interception from the canopy, the total incoming radiation can be expressed as:

${\displaystyle (1)\,\,\,\,={\underset {az}{\iint }}(a,z)\,da\,dz}$

where (a, z) is amount of radiation (direct plus diffuse) incoming from the direction (a, z) over the simulation time step. In the presence of a canopy, the incoming radiation is attenuated as:

${\displaystyle (2)\,\,\,\,={\underset {az}{\iint }}(a,z)\cdot e^{-k\cdot LA(a,z)}\,da\,dz}$

where LA is the leaf area intercepted by the ray of direction (a, z) and k is a parameter, for which an average value of 0.5 has been assumed. By subdividing the sky hemisphere into a finite number n of solid angles, each associated with azimuth a and zenith z of its centroid, <IATot> can be estimated by Monte Carlo integration as:

$I}_{ATot}> = \frac{1}{n} \sum_{i=1}^n \frac{f([a,z]_i)}{ p([a, z]_i)} \\ \\ \\ \mathrm{where } \quad f([a,z]_i) = ([a,z]_i)\cdot e^{-k \cdot LA([a, z])_i)$

Using importance sampling, p([a, z]i) should follow as closely as possible the expected form of f([a, z]i) to reduce the error of the estimate for the same value of n. The value of <I>([a, z]i) in equation 3 is calculated for each time step using the approach of Brunner (1998): diffuse radiation is expressed as function of z under standard overcast sky (SOC) conditions, while direct radiation is expressed as a function of a and z from the corrected number of sun positions. The equations of Brunner (1998) permit the calculation of the spatial distribution of a normalised above canopy light (ACLa,z), but since by definition

${\displaystyle (4)\,\,\,\,ACL([a,z]_{i})={\frac {([a,z]_{i})}{}}}$

then it is possible to calculate

${\displaystyle (5)\,\,\,\,<{\hat {I}}_{ATot}>=\cdot {\frac {1}{n}}\sum _{i=1}^{n}{\frac {ACL([a,z]_{i})\cdot e^{-k\cdot LA([a,z]_{i})}}{p([a,z]_{i})}}}$

where <ITot> is obtained either by direct measurements from comparable unobstructed sites, from existing data sets (Pinker & Laszlo 2000;Scharmer & Greif 2000), or from remote sensing based estimates (Frouin & Pinker 1995). The next step in the model is to define a distribution of p([a, z]i) which ensures the largest reduction of the error in the estimate for the same value of n. To this end, the distribution of p should match the distribution of ACL as well as the expected spatial distribution of LA([a, z]i). We propose the form:

${\displaystyle (6)\,\,\,\,p([a,z]_{i})={\bigg (}{\frac {w1}{w1+w2}}\cdot {\frac {p1(ACL([a,z]_{i}))}{\sum _{i=1}^{n}p1(ACL_{i})}}+{\frac {w2}{w1+w2}}\cdot {\frac {p2(eLA([a,z]_{i}))}{\sum _{i=1}^{n}p2(eLA_{i})}}{\bigg )}}$

where p1 and p2 are the probabilities associated with the two determinants of incoming light, w1 and w2 are weights that can be used to modulate the relative importance of ACL and LA, and eLA is a function that expresses the expectation of LA for a given a and z. The distribution of ACL over the sky hemisphere, calculated using the approach of Brunner (1998) prior to dynamic stand modelling, is the same for every point in the simulated stand, so ACL can be used directly by setting:

${\displaystyle (7)\,\,\,\,p1(ACL([a,z]_{i}))=ACL([a,z]_{i})}$

The value of p2 should match the expected variability in LA, increasing its value where LA is more likely to be variable. It will be modelled as a function of z, since there are no a priori reasons to expect a change in variability of LA with a. Three different probability distributions will be employed:

1. A function for regeneration plots, where p2 declines with increasing z: the sampling is concentrated in the upper 45 degrees of the hemisphere, since in gaps and small openings lateral light is much more likely to be intercepted by the surrounding overstorey
2. A function for saplings and intermediate trees which weighs equally lateral rays and rays coming from above, to assess both lateral competition and the effect of the canopy of dominant trees
3. A function for dominant and co-dominant trees, where sampling is concentrated in the lateral portion of the hemisphere, since the sky above is likely to be unobstructed.

The definitive format of the sampling probability functions, as well as the sampling intensity, will be determined through sensitivity analysis, considering trade-offs between precision and computational efficiency. Known the sampling intensity and the sampling probabilities, a set of rays are randomly selected for each simulation time step, and used in the ray-tracing routine to calculate the light environment at the top of every regeneration plot as well as each individual tree. Figure 1 shows the relationship between the part of the algorithm that has to be calculated only once per time step and the parts that must be calculated for every individual point. The ray-tracing algorithm uses the position of the point and the stand map to calculate the number of crowns intercepted by each sampling ray. It is assumed that

${\displaystyle (8)\,\,\,\,LA([a,z]_{i})=\sum _{j}(l_{j}\cdot hits([a,z]_{i})_{j})}$

where hitsj is the number of crowns of the species j intercepted by the ray and lj is a species-specific leaf area coefficient (Canham et al. 1994;Canham et al. 1999). This is equivalent to assuming that the foliage of a tree is concentrated in a layer of constant thickness and uniform density on the surface of the crown. The light environment at point (x, y, z) is then calculated by substituting the result of Equation 8 in Equation 5.