Consideration of radiation absorption by stems in forests for microclimate modeling

Forest canopy models are used to simulate biosphere–atmosphere coupling in global climate models, as well as the microclimate influences of forests at the stand scale. It has recently been shown that wood structures store significant heat following radiation absorption and impact air temperature diurnal patterns inside canopies. Yet, radiation absorption by woody stems is not fully considered in current models. Here we modify the radiative transfer component of the CanVeg2 multilayer canopy model to include radiation absorption by woody stems. We evaluate the model modifications by comparing estimates against a 3D ray tracing radiative transfer model parametrized using ground lidar measurements, and against tower observations of albedo in four broadleaf forests. We found a very good agreement between the 1D and 3D models, and a good agreement between models and observations. Our approach provides a tractable and computationally efficient implementation of radiation absorption by woody stems to calculate biomass heat storage in canopy models.

1 Introduction

The microclimate of forests in terms of air temperature and humidity often varies vertically from the soil surface to the canopy top, and these conditions can significantly differ from those above the canopy (Bonan, 2023). Such buffering effects have implications for organisms living within the canopy and near the soil surface (Noble et al., 2026). The scalar variations within a canopy are governed by an energy balance between radiation, conduction and latent and sensible heat fluxes which occur at the surface of – and within – leaves, stems, and soil. Woody stems typically absorb a significant portion of incoming near infrared radiation in dense forests (Béland, 2025), and most of the absorbed radiation is transferred as heat within stems during the day and is released at night. For the purpose of modeling forest microclimate, it is thus relevant to consider stems in the energy balance, and scalar vertical profiles through a multilayer canopy model (Bonan et al., 2021). Further, the radiative forcing on stems leads to potentially significant amounts of heat storage which may play a role in the energy budget closure of eddy covariance flux tower measurements (Leuning et al., 2012).

Radiation absorption in forest canopies is commonly modeled using one dimensional radiative transfer models like two-stream (Dickinson, 1983; Sellers, 1985) and the Norman model (Norman, 1979). They compute the portion of solar energy that is reflected (albedo) and allocate the absorbed energy portion to leaves (often distinguishing sunlit and shaded) and soil. The model computations are based on canopy structure information like the leaf area index (LAI) and leaf angle distribution functions, as well as leaf level optical properties (typically separated between visible (PAR) and near-infrared (NIR) radiation). Neither model currently accounts for a heterogenous mixture of absorbing and scattering elements (e.g., leaves and wood) in the canopy. At global scales, the two-stream implementation in the Community Land Model (CLM) uses plant area index (leaf area plus stem area) to calculate radiative transfer based on a weighted combination of leaf and stem optical properties. The separate absorption of radiation by leaves and stems is not considered. In Béland et al. (2026), we circumvented this problem by using a 3D radiative transfer model, but this is not practical for running microclimate simulations using models like CanVeg2 at sites where 3D canopy structure description is not available. Hence, we seek to develop a more generalizable approach to modeling radiation absorption by stems for use in microclimate models.

The absorption of radiation by stems is challenging to model, mainly because they are highly heterogeneous in space, both in terms of size, orientation and density, and because they are exposed to different levels of solar radiation depending on their position within the canopy. The advent ground lidar provides a means of mapping the area density of leaves and wood separately in three-dimensional space (Béland and Kobayashi, 2024; Béland, 2025). This enables the estimation, through ray tracing radiative transfer modeling, of the radiative forcing on stems of different sizes located at different heights above ground. Ground lidar surveys and data processing require specialised equipment and expertise which limits wide applicability, but provide a valuable means of validating the estimation of radiation absorption by stems from a 1D radiative transfer model and 1D leaf and wood area index profiles.

Here, we present modifications to the radiative transfer component of the CanVeg2 multilayer canopy model (Béland et al., 2026) to account for the absorption of radiation by wood structures. CanVeg2, like the original CanVeg, can use the Norman (1979) radiative transfer. We used canopy structure descriptions derived from ground lidar and a 3D ray tracing radiative transfer at four broadleaved deciduous sites to validate the modifications. The outputs from the 3D ray tracing and Norman models provide a benchmark for other canopy models.

2 Models description

2.1 Modified Norman's model and its solution using a tridiagonal system of equations

The Norman model consists in a system of equations describing the scattering and absorption of direct and diffuse radiation by canopy layers each having a leaf area index (ΔLAI). The model does not consider multiple scattering within a layer, so thin layers should be used. The model equations mainly describe the fluxes of transmitted and absorbed direct and diffuse radiation by each layer as the radiation travels downward and upward (following scattering). Both the transmission and absorption are calculated following the Beer-Lambert exponential law involving the G function and foliage clumping which respectively describe the leaf inclination angle distribution and the aggregation of foliage at different spatial scales (Monsi and Saeki, 1953; Ross, 1981; Nilson, 1971). Béland and Baldocchi (2021) presented the model equations modified to consider foliage clumping in the canopy layers.

The modifications to the Norman model presented here essentially consist in combining the projected area of leaves and stems in the Beer-Lambert equation, and reworking the algebra involved in the tridiagonal system of equations to include the stem contribution to transmission and absorption fluxes. The modifications for the transmission portion are straightforward; the projected surface areas for leaves and stems can be summed in the Beer-Lambert equation (see Supplement). A layer's intercepted fluxes are assigned to leaves and stems by calculating their respective fractions of projected area within a given layer. In other words, the contribution of one element (leaves or wood) to the interception of radiation from a given direction in a given layer is proportional to the ratio of the projected area of the element in this direction to the sum of the projected areas of both elements. Mathematically, for shortwave direct beam radiation, this is expressed for leaves as:

$$\begin{array}{}\text{(1)}& {\mathrm{frac}}_{\mathrm{b},\mathrm{l},i}={\displaystyle \frac{G_{\mathrm{l},i}{\mathrm{\Omega }}_{\mathrm{l},i}\mathrm{\Delta }{\mathrm{LAI}}_i}{G_{\mathrm{l},i}{\mathrm{\Omega }}_{\mathrm{l},i}\mathrm{\Delta }{\mathrm{LAI}}_i+\mathrm{2}\cdot G_{\mathrm{w},i}{\mathrm{\Omega }}_{\mathrm{w},i}\mathrm{\Delta }{\mathrm{WAI}}_i}}\end{array}$$

And for wood as:

$$\begin{array}{}\text{(2)}& {\mathrm{frac}}_{\mathrm{b},\mathrm{w},i}={\displaystyle \frac{\mathrm{2}\cdot G_{\mathrm{w},i}{\mathrm{\Omega }}_{\mathrm{w},i}\mathrm{\Delta }{\mathrm{WAI}}_i}{G_{\mathrm{l},i}{\mathrm{\Omega }}_{\mathrm{l},i}\mathrm{\Delta }{\mathrm{LAI}}_i+\mathrm{2}\cdot G_{\mathrm{w},i}{\mathrm{\Omega }}_{\mathrm{w},i}\mathrm{\Delta }{\mathrm{WAI}}_i}}\end{array}$$

where frac${}_{\mathrm{b},\mathrm{l},i}$ and frac${}_{\mathrm{b},\mathrm{w},i}$ are the fraction of intercepted beam radiation by leaves in layer i, and by wood respectively. ΔLAI_i and ΔWAI_i are the leaf area index and wood silhouette area index in layer i. Ω_l,i and Ω_w,i are the leaf and wood clumping indices in layer i. And G_l,i and G_w,i are the leaf and wood projection factors for layer i. The same principle is applied to the intercepted diffuse radiation. The fractions of absorbed and scattered radiation by leaves and wood from an intercepted flux in a given layer are then determined by multiplying the intercepted flux by frac${}_{\mathrm{b}/\mathrm{d},\mathrm{w}/\mathrm{l},i}$ and the optical properties of the leaves/wood. The full set of equations adapted from the numerical methods presented in Bonan (2019) to solve Norman's model equations are provided in the Supplement for shortwave and longwave radiation.

2.2 The clumping of stems and branches

The foliage clumping index is used in radiative transfer models to account for the effect of leaf spatial distribution on radiation interception (Chen et al., 2012). It is reasonable to suggest that stems and branches are also clumped and it is useful to describe this clumping for radiative transfer purposes. Kucharik et al. (1998) suggested a non-randomness factor for branch location in canopy space. Since the concept of foliage clumping is linked to the non-random distribution of elements in space through the Markov model, it is appropriate to apply the concept to stems and branches.

To estimate a stem clumping index, we used light transmission data collected from a tram system in a temperate deciduous forest near Oak Ridge, Tennessee, in the 1980s (Baldocchi et al., 1985). Some of the data was collected in leaf-off conditions on 20, 25 and 26 January 1981 from downward PAR sensors traveling back and forth along the 30 m long tram at about 2 m above ground and above the canopy. The data thus represents a horizontally averaged measure of transmitted PAR which can be used to invert Beer's law and estimate the clumping of wood structures:

$$\begin{array}{}\text{(3)}& \mathrm{\Omega }={\displaystyle \frac{-\ln \left(P_{\mathrm{0}}\right)\cdot \cos \left(\mathit{\theta }\right)}{\mathrm{WAI}\cdot \mathrm{2}\cdot G\left(\mathit{\theta }\right)}}\end{array}$$

where P₀ is the fraction of transmitted light, θ is the sun zenith angle, G(θ) is the element projection factor in the direction θ, and WAI is the wood silhouette area index. A value of 0.5 for G was used, which is reasonable for wood structures as reported in Béland (2025). The factor 2 appearing in the denominator reflects the non-flat nature of wood structures and the use of the wood silhouette area index to define wood structures, as discussed in Béland (2025). The value for WAI of 0.6 for this site reported in Hutchison et al. (1986) was used. One needs to consider, however, that the near ground level measurements include some scattered light, and since the wood bark has a higher scattering coefficient than leaves, the contribution of this scattering to the PAR received at the ground level is expected to increase light transmission in a non-negligible way.

Applying Eq. (3) to the three days of light transmission data between noon and 02:00 p.m. when the sun zenith angle was between 30 and 35° resulted in an average stem clumping index of about 0.65. We estimated using 3D ray tracing modeling that the effect of the increased light transmission caused by light scattering leads to a decrease of the stem clumping index by about 0.1. Hence, we estimate of the wood clumping index in a temperate deciduous forest to be 0.75; this value was used here for Ω_w without any vertical variability.

2.3 The clumping of leaves

Foliage clumping – first defined through the effect of leaf aggregation level on light transmission by Nilson (1971) – has long been known to lower canopy albedo (Dickinson, 1983; Ni-Meister et al., 2010), but one dimensional radiative transfer models have yet to fully consider this effect. It can be considered through the probability of recollision between a photon and canopy parts through its effect on the scattering coefficient, proposed by Ross (1981) and integrated into more complex modeling schemes using tree shoot measurements by Smolander and Stenberg (2003). The relationship between leaf clumping index and recollision probability has been studied by Stenberg et al. (1994) and Smolander and Stenberg (2003), who found that the recollision probability is equal to one minus the clumping index at the shoot scale. Their work focused on needleleaf shoots, but it is reasonable to suggest the theory applies to clumping in broadleaf shoots as well.

The relationship between shoot clumping index, or shoot shading factor (equal to four times the silhouette to total area ratio ($\mathrm{4}\stackrel{\mathrm{‾}}{\mathrm{STAR}}$) in Stenberg, 1996 and Smolander and Stenberg, 2003), and recollision probability is an approximation used to derive the scattering properties of a shoot as a function of its structure. The 3D radiative transfer model used here (FLiESvox, presented in Sect. 3.3) uses this approximation to determine the occurrence of multiple scattering within the elementary volume (voxel cube), and we maintained this relation in the Norman model to consider multiple scattering within layers. This was done by modifying a layer's scattering coefficient following the approximation proposed in Smolander and Stenberg (2003) (see details in Supplement).

It remains unclear how shoot leaf area density influences the relation between shoot shading factor and recollision probability. We choose to use this theory here to maintain consistency between the 3D model and the modified 1D Norman model, and because we believe it is useful in improving the correspondence between modeled and measured albedos, as well as the accuracy of the amount of radiation absorbed by stems, which has consequences for stem heat storage and the canopy energy budget. We do note, however, that the physical basis and parametrization of the theory remain uncertain.

3 Experiment Description

We used canopy structure data derived from ground lidar at four sites along a canopy structure gradient, and 3D ray tracing radiative transfer modeling to assess the correspondence with the modified Norman radiative transfer model. We also use field measured leaf optical properties combined with tower measured albedo to assess correspondence with the radiative transfer modeling estimates of canopy albedo.

3.1 Sites Description

Four sites that varied in vertical structure following a gradient – from bottom-heavy to top-heavy canopies – were used. Each site was surveyed using a ground lidar Riegl VZ-400 from different positions forming a very dense grid of 121 scan positions to minimize occlusions and provide a complete representation of the leaves and wood structures within 60 m × 60 m plots. The lidar data was processed using methods presented in Béland and Kobayashi (2021) to yield a 3D voxel grid of leaf area density, each voxel having 30 cm in side length. The wood structures were also modeled using a voxel grid from similar methods presented in Béland (2025) to yield wood silhouette area density.

The first site is at the Harvard Forest Environmental Monitoring Site (EMS), located in Petersham, Massachusetts, USA (42°32^′ N, 72°10^′ W, elevation 340 m). The mean canopy height is about 25 m with a leaf area index (LAI) estimated at 4.8 and dominated by red oak (Quercus rubra) and red maple (Acer rubrum). The leaf angle distribution is planophile and constant vertically (Béland and Kobayashi, 2021). Leaf spectral properties measurements collected by Bartlett et al. (2011) using an ASD FieldSpec 3 (Analytical Spectral Devices Inc., Boulder, Colo.) connected to an integrating sphere (SphereOptics, Concord, N.H.) were available for this site. The selected plot is located between two eddy covariance flux towers, one operated by Harvard University and the other by the National Ecological Observatory Network (NEON).

The second site is in the Morgan-Monroe State Forest located near Bloomington, Indiana, USA (39°19^′ N, 86°25^′ W, elevation 250 m). The mean canopy height is about 30 m with a LAI estimated at 5.1 and dominated by sugar maple (Acer saccharum), white oak (Quercus alba), red oak (Quercus rubra) and tulip poplar (Liriodendron tulipifera). The leaf angle distribution is planophile up to about 30 m and then transitions to uniform in the emergent crowns (Béland and Kobayashi, 2021). This site has an eddy covariance flux tower and is part of the Ameriflux network.

The third site is at the Smithsonian Environmental Research Center (SERC) located near Edgewater, Maryland, USA (38°53^′ N, 76°33^′ W, near sea level). The mean canopy height is 33 m with an LAI of 6.2 and dominated by American beech (Fagus grandifolia) and tulip poplar (Liriodendron tulipifera). The leaf angle distribution is planophile up to about 30 m and transitions to erectophile in the emergent crowns (G. Parker, personal communication, 2019). Leaf optical properties measured by NEON using a ASD field spectrometer connected to a leaf clip (Analytical Spectral Devices Inc., Boulder, Colo.) in 2016 were available (NEON, 2022). This site also has an eddy covariance flux tower operated by NEON.

The fourth site is in the Pasoh Forest Reserve of the Forest Research Institute Malaysia located near Simpang Pertang, Malaysia (2°58^′ N, 102°18^′ E, elevation 140 m). The mean canopy height is 35 m with an LAI of 7.3 (Kira et al., 2013) and composed of mixed dipterocarp trees. The leaf angle distribution data for this site was observed from 30 m above ground and fits the uniform distribution (Y. Kosugi and S. Takanashi, personal communication, 2019), it is likely that leaves tend towards a planophile distribution below 30 m and erectophile above. This site has an eddy covariance flux tower and is part of FLUXNET (Pastorello et al., 2020).

3.2 Lidar data and processing

The ground lidar collected at each site allows distinguishing between signal received from leaves and wood, which is a very useful feature of ground lidar. Methods to achieve this classification are presented in Béland and Baldocchi (2021); they rely both on the intensity of the reflected laser light and on the organization of points in 3D space. The different canopy structures are represented using voxel arrays where each voxel contains the voxel's leaf area density, leaf angle distribution function (which can change vertically as described in the previous section), and wood silhouette area density. This information can in turn be used to generate inputs for 1D or 3D radiative transfer models (for 1D models, vertical profiles of leaf and wood silhouette area density are generated by horizontally integrating the leaf and wood area of all voxels in a given layer).

The 3D voxel arrays were generated using methods presented in Béland et al. (2014), Béland and Kobayashi (2021) and Béland (2025). The leaf area density within each voxel is calculated using a ray tracing algorithm called VoxLAD_RT (Béland, 2026a), which traces the path of each laser pulse emitted by the lidar instrument as they travel through and are intercepted within a given voxel. The contact frequency method developed by Warren Wilson (1960) is used to estimate the leaf area density and wood area density from information on the leaf/wood orientation, the number of laser pulse intercepted and the length of their travel paths within the voxel. The within voxel foliage clumping profiles were established from lidar and direct measurements made on tree branches collected from different heights above ground. Each branch is placed within a cubic frame with open sides and scanned using a ground lidar instrument, each leaf located within the cubic frame is then collected to calculate de leaf area density of the sample. The lidar measurements provide the leaf orientation and light transmission through the cubic volume, which enable calculating the clumping factor (Béland and Baldocchi, 2021).

3.3 3D ray tracing radiative transfer modeling

The 3D ray tracing radiative transfer model FLiESvox (Kobayashi and Iwabuchi, 2008; Béland and Kobayashi, 2024) was used here to estimate the amount of radiation in the PAR and NIR absorbed in the different vertical canopy layers (30 cm in thickness) by leaves and wood for comparison against the modified Norman model (a simulation example is shown in Fig. 1). FLiESvox can be considered a hybrid model in the sense it uses the voxel as the basic structure element within which leaf clumping is determined by field observations made at the voxel scale (see Béland and Baldocchi, 2020, 2021), while larger scale leaf clumping across the site domain is determined by the spatial distribution of the voxels and their leaf area densities. In the case of wood, the wood area is randomly distributed inside a voxel, and the clumping of wood structures arises only from the distribution of voxels in 3D space. Any voxel can contain both leaf and wood area, and the probability of a photon entering the voxel to interact with a leaf or wood is partly determined by the leaf area density, the within voxel leaf clumping, the wood area density, and the angular distribution function of both elements inside the given voxel. Upon an interaction with clumped leaves, if the photon is scattered (which is probabilistically determined by the leaf optical properties), the recollision probability theory of Smolander and Stenberg (2003) is used to calculate the probability that the photon will interact with another leaf inside the voxel as one minus the leaf clumping factor associated with the voxel's height above ground (derived from the shoot level clumping vertical profile mentioned in Sect. 2.2). In this 3D model, the recollision probability is used to calculate the probability a scattered photon will intercept another leaf within the voxel because leaves are not explicitly positioned in space within voxels. At scales larger than the voxel, the ray tracing algorithm calculates the probability of photon interception based on leaf area density, optical properties, orientation, and the photon travel path length within the voxel.

Figure 1 Examples of a 3D voxel-based mapping of absorbed radiation (in W m⁻²) by leaves (green) and wood (brown) in PAR (left) and NIR (right) as simulated by the FLiESvox 3D raytracing radiative transfer model at the Harvard Forest site.

3.4 Radiative transfer models parameters

Two sets of model runs were performed towards two objectives: (1) compare 1D and 3D radiative transfer simulations to confirm the effect of adding new functionalities like increased multiple scattering from leaf clumping or radiation absorption by wood, and (2) compare tower observation of albedo in different illumination conditions to radiative transfer model estimates.

In model runs aimed at comparing simulated and tower-based observation of albedos, the observed incoming radiation diffuse fractions in PAR and NIR as well as the sun zenith angle at the time of observations (every 1 h at EMS and Morgan Monroe and 30 min at SERC and Pasoh) were used in the Norman model. Since it would be computationally impractical to run the 3D FLiESvox model of each illumination conditions in the long time series used, we used a look up table containing FLiESvox radiative transfer simulation results for 9 different sun zenith angles and 9 different incoming radiation diffuse fractions, for a total of 81 PAR model runs and 81 NIR runs to cover different illumination conditions. The nearest values in the look up table to the actual diffuse radiation fractions and sun angles were used.

In model runs aimed at comparing the 1D and 3D simulations, both models were reconfigured to evaluate the effect of the different modifications made to the Norman model. The Norman model configurations were: (1) not considering foliage clumping nor wood structures, (2) considering foliage clumping only in terms of its effect on light penetration (recollision probability is not considered) and wood structures were excluded (this corresponds to the version of Norman's model presented in Béland and Baldocchi, 2021), (3) same as (2) with the addition of wood structures, and (4) all modifications considered: foliage clumping effect on light penetration and recollision probability, and wood structures. All these model runs used a fixed illumination configuration of 35° sun zenith angle, 20 % diffuse incoming radiation, and a total incoming radiation flux of 500 W m⁻². In order to evaluate the fit between the 1D and 3D models before the modifications to the Norman model presented here were applied, the FLiESvox model was reconfigured to exclude the wood structures and the recollision probability, hence considering only the effect of foliage clumping on light penetration.

Both the modified Noman model and the 3D ray tracing model were parametrized using the same values for leaf optical properties and illumination conditions. A leaf reflectance value in the PAR of 6 % was used at all sites, this value corresponded to field measurements at Harvard Forest and SERC made by Bartlett et al. (2011) and Neon (2022), and is 2 % lower than the average value for broadleaf deciduous trees species published in Majasalmi and Bright (2019). A leaf transmittance value in the PAR of 7 % (corresponding to the value measured at the Harvard Forest site) was used at all sites except SERC, where a value of 5 % corresponding to the measurement of NEON (2022) (Majasalmi and Bright (2019) reported an average value of 6 %). In the NIR leaf reflectance and transmittance were set at 45 % each at all sites, corresponding to the measurements from NEON (2022). This yields a 10 % absorbance, which is 5 % lower than the average value of broadleaf deciduous tree species published in Majasalmi and Bright (2019), but their database did not include any of the species present at the sites used here. We did not consider the NIR data collected from Bartlett et al. (2011) since for several samples the sum of reflectance and transmittance was higher than 100 %, suggesting potential issues with field data collection methodology. The stem reflectance in the PAR was set to 21 % and in the NIR at 49 %, corresponding to the values suggested by Majasalmi and Bright (2019) for broadleaf deciduous trees. The soil reflectance was set to 12.2 % in the PAR, and 21.4 % in the NIR, corresponding to average values used in CLM (Bonan et al., 2018).

The leaf area density profiles used in the Norman model came from summing the leaf area density values for all voxels in each vertical layer, and the vertical resolution for the profile corresponded to the voxel height of 30 cm. The vertical profiles of foliage clumping used in the modified Norman model were those published in Béland and Baldocchi (2021) at a vertical resolution of 1 m.

3.5 Tower based observations

The tower-based albedo observation were made using cosine receptors, and the periods used were for the month of July between 2017 and 2020 for Harvard Forest (NEON, 2021a), 2015–2020 for Morgan Monroe (Novick and Phillips, 2023), 2018–2021 for SERC (NEON, 2021b), and for the month of February between 2003 and 2009 for Pasoh (Kosugi and Takanashi, 2009). The diurnal patterns of albedo were calculated as the average albedo for each hour (or half-hour at the SERC and Pasoh sites) over the month. The incoming and outgoing shortwave measurements were available at all sites, and the PAR measurements were available at the Harvard Forest and SERC sites, allowing the calculation of NIR albedo (calculated from the shortwave and PAR fluxes).

4 Results

4.1 Comparison between models

A first comparison of reflected and absorbed flux estimates between different versions of the Norman model against the 3D ray tracing model FLiESvox is shown in Fig. 2 for the Harvard Forest site (data for other sites are provided in Figs. S1 and S2 in the Supplement and follow similar trends). Shown are the original Norman model from Norman (1979), the Norman model with the addition of foliage clumping as presented in Béland and Baldocchi (2021), the modifications to the Norman model presented here to include the interception of radiation by wood structures, and the introduction of multiple scattering within layers of clumped leaves.

These results show that in the NIR (top panel), the original Norman model (no clumping, no wood, no recollision) has very reasonable proportions of radiation absorbed by shaded and sunlit leaves, as well as by the soil compared to the 3D ray tracing model, but has a much higher albedo. The addition of foliage clumping in the Norman model decreases albedo and increases the amount of radiation reaching the soil – since light is able to penetrate deeper in the canopy –, while proportions of light absorbed by leaves are reduced. Adding the wood structures to the Norman model further reduces albedo and soil absorption, and slightly reduces leaf absorption. The addition of multiple scattering in clumped layers significantly reduces albedo and increases leaf absorption. The overall agreement between the 3D ray tracing estimates and the modified Norman model which includes wood structures and multiple scattering is considered high and is further analyzed in Fig. 3.

In the PAR wave band, the original Norman model results in a high proportion of light absorbed by sunlit leaves, and little light reaches the soil. Adding foliage clumping to the model decreases the amount of radiation absorbed by sunlit leaves and increases the radiation absorbed by shaded leaves, establishing proportions very close to the 3D ray tracing model values. The addition of wood structures further decreases radiation absorption by sunlit leaves and soil. Adding multiple scattering has relatively little effect since leaf-level scattering is very low in the PAR.

Figure 2 Budget for the proportions of NIR (top) and PAR (bottom) reflected, and absorbed by sunlit and shaded leaves, by wood and soil, simulated using the Norman model with four different configurations and the FLiESvox 3D ray tracing model. All simulations are for sun zenith angle of 35°, diffuse fraction of 20 % and incoming radiation of 500 W m⁻². Shown here are results for the Harvard Forest site (equivalent figures for all other sites are provided in the Supplement).

Figure 3 Correspondence between the total radiation fluxes reflected and absorbed for the modified Norman and FLiESvox models in the NIR (top) and PAR (bottom) for all sites. The Morgan Monroe site is in orange, Harvard Forest in yellow, Pasoh in green, and SERC in blue. The red line is 1:1.

Figure 3 provides a direct comparison in proportions of absorbed radiation when using the modified Norman vs. the 3D ray tracing model across all four sites used. In the NIR the agreement is good with no apparent bias. In the PAR, three of the four sites show an overestimation of radiation absorbed by wood.

Previous results concerned vertically integrated absorbed radiation fluxes, Fig. 4 shows the vertically resolved profiles of absorbed radiation for the Harvard Forest site. This site is perhaps where the correspondence between the 1D and 3D models is best, data for other sites is provided in the Supplement (Figs. S3 to S8) and the correspondence is generally considered very good across sites, with a mean RMSE per vertical leaf layer of 0.59 W m⁻² in PAR and NIR. The main source of differences is believed to originate from the way foliage clumping is considered in the Norman model. This is based on near perfect agreement between 1D and 3D models when considering only a turbid medium without any clumping (data not shown), while the vertical profiles from model runs considering only foliage clumping (no wood structures or recollision probability) show a divergence with a mean RMSE per leaf layer of 0.72 W m⁻² (data shown in the Supplement, Figs. S9–S12). This may be explained by some 3D effects being lost in 1D model, and the fact that even with recollision probability turned off in the 3D model, multiple scattering inside a layer remains possible because of the horizontally resolved heterogeneity, which is not the case in the Norman model. This leads to increased absorption by leaves especially in the NIR.

Figure 4 Vertical profiles of absorbed radiation by sunlit leaves (top), shaded leaves (middle), and wood (bottom) in the PAR (left) and NIR (right) for Norman (blue lines) and FLiESvox (red lines) models. Data shown here is for the Harvard Forest site, data for all sites is provided in the Supplement.

4.2 Comparison to tower albedos observations

Over the entire shortwave, the 3D model estimates were slightly lower than the 1D model estimates at all sites as shown in Fig. 5. The tower albedo observations were lower than the modeled values at all sites. The shapes of the diurnal patterns were similar between modeled values and observations. The results show an improvement in the correspondence between modeled and measured albedo in the modified vs original Norman model, which is attributed to radiation absorption by stems and the consideration of multiple scattering in canopy layers.

Figure 5 Tower observed and modeled using the modified Norman (in yellow), Norman model version presented in Béland and Baldocchi (2021) (here labeled Original Norman, in red), and FLiESvox (in green) models for shortwave albedo at all sites. Diurnal patterns are averages for each hour at Harvard Forest and Morgan Monroe, and half hour at SERC and Pasoh, for the month of July at all sites except Pasoh which is for February.

To characterize the sensitivity of the modeled albedos to leaf optical properties, we performed an additional run of Norman's model for the Morgan Monroe site using a leaf reflectance of 40 % instead of 45 % in the NIR, while maintaining leaf transmittance to 45 %, reducing leaf absorptance from 10 % to 15 %. The result, shown in Fig. 6, indicates a decrease of broadband albedo by about 4 %, which attests to the impact of leaf optical properties on albedo, especially in the NIR and when leaves are clumped at the canopy top which leads to increased multiple scattering of photons.

Figure 6 Demonstrating the effect of a change in leaf optical properties on simulated shortwave albedo for the Morgan Monroe site. In yellow is shown the Norman model simulation output using the same value for leaf reflectance in the NIR 45 % as elsewhere in this study, and in green using leaf reflectance of 40 %.

5 Discussion

5.1 Comparison between models

Both the 1D and 3D models make use of common principles in the treatment of multiple scattering within a volume of clumped foliage, and they both used the same leaf area density vertical profiles and similar foliage clumping profiles, it is thus expected that they have a fair level of agreement. They also have several independent features, like the nature of their absorbed radiation calculations (ray tracing vs. equations with a numerical solution) and the full consideration of the 3D placement of leaves and wood in the 3D model. Given the importance of independent elements, the fit between the 1D and 3D models is considered very good, agreeing to about 12 W m⁻² in the PAR in terms of RMSE, and about 8 W m⁻² in the NIR (see Fig. 3) over all vertically integrated layers of absorbed fluxes at all sites. Vertical profiles of absorbed fluxes are also in very good agreement, the average RMSE across sites for all vertical layers of PAR absorbed by sunlit leaves was 0.9 W m⁻², 0.1 W m⁻² for PAR absorbed by shaded leaves, 0.2 W m⁻² for NIR absorbed by sunlit leaves, 0.1 W m⁻² for NIR absorbed by shaded leaves, 0.2 W m⁻² for PAR absorbed by wood and 0.2 W m⁻² for NIR absorbed by wood (see Figs. S3 to S8 in the Supplement).

How well the modified Norman model would perform in open canopies like savannas and in needleleaf forests remains to be assessed. We did however perform an additional model run to assess the correspondence between 1D and 3D models in the leaf-off case, i.e. where leaves are absent and only wood structures remain. Both models strongly agreed (see Figs. S13 and S14), suggesting the modified Norman model presented here is also valid for simulating radiation absorption by wood during periods of the year when leaves are absent. To assess the modified Norman model sensitivity to the wood clumping factor and diffuse radiation fraction, we performed additional model runs where the clumping factor is decreased from 0.75 to 0.6 and increased to 0.9, and the diffuse fraction is decreased from 20 % to 50 % and 90 %. The results shown in Figs. S15 to S16 in the Supplement indicate that decreasing the wood clumping factor to 0.6 (thus increasing the level of wood aggregation) decreased the radiation flux absorbed by stems by about 16 % in a consistent manner across sites and wavebands. Increasing the wood clumping factor to 0.9 also increased the radiation flux absorbed by stems by about 16 % across sites and wavebands. Decreasing or increasing the wood clumping factor had a very modest effect on simulated canopy albedos, with average differences around ±0.5 % across sites in the NIR, and ±0.03 % in the PAR. The value of 0.75 for the wood clumping factor used here was derived from one site, and it may be relevant considering its noticeable impact on stem radiation absorption to pursue further research on estimating this factor in different forest types. Results also show that diffuse radiation fraction did not impact noticeably the amount of radiation absorbed by stems in neither PAR or NIR.

The FLiESvox model is an offspring of the FLiES model which has been extensively tested in the RAMI radiative transfer benchmarking exercise (Widlowski et al., 2015, 2013; Lanconelli et al., 2025). FLiESvox simulations have been compared against airborne hyperspectral observations (Béland and Kobayashi, 2024). These evaluation results lead us to have high confidence in this model as a gold standard for benchmarking simpler models. All input data and model outputs for the test cases presented here are made available online (see data availability section) for future model modification evaluation exercises.

5.2 Comparison to tower albedo observations

Measurements of absorbed radiation by wood, sunlit and shaded leaves are not feasible at the plot scale, and measurements of absorbed radiation by soils require an approach to account for the horizontal heterogeneity of radiation fluxes reaching the ground (sometimes done using a tram system (Hutchison et al., 1986)) and are rather rare. Hence, of all components simulated by the models, canopy albedo is the only one for which measurements are commonly available, and consequently perhaps the only measurement that can be used to directly evaluate the vertically integrated radiation fluxes estimated from a canopy radiative transfer model.

Figure 5 shows the correspondence between model estimates and tower-based broadband albedos (including both PAR and NIR). Similarly to the difference between model estimates and NIR albedo, both the 1D and 3D models overestimate the tower-based broadband albedo by between 4 % and 7 %, with larger differences at Morgan Monroe and Pasoh. Past study by Baldocchi et al. (1985) have shown poor agreement between Norman's model and scattered radiation. They noted that foliage clumping had a role in the overestimation of modeled values compared to observations from a flux tower. Here, even though the models overestimate albedo, the addition of radiation absorption by wood and multiple scattering in clumped volumes using the recollision probability theory considerably improved the match to ground observations compared to the Norman model version presented in Béland and Baldocchi (2021) by lowering NIR albedo (see also Fig. 2) while retaining reasonable correspondence in leaf level optical properties with regards to field measurements. Further, the “U” shape of the diurnal albedo pattern (higher albedos at lower sun elevation angles) was successfully reproduced by the modified Norman model.

In this study, we assumed that the recollision probability theory at the shoot scale, initially developed for needleleaf shoots (Smolander and Stenberg, 2003), also applies to broadleaf shoots, and that the relation between shoot clumping index and recollision probability is a good enough approximation to account for multiple scattering in clumped element volumes. This is based on recent work showing that broad leaves can be clumped at the shoot scale (Béland and Baldocchi, 2021), and work by Stenberg et al. (1994) and Smolander and Stenberg (2003) who relate the clumping index to recollision probability. In the light of the fit between modeled and observed albedos obtained here, it appears that further experiments are needed to test this relationship specifically for broad-leaf shoots of different morphology and leaf area density. This future research could make use of ray tracing models where leaves are individually explicitly represented using meshes, like in the DART model (Gastellu-Etchegorry et al., 2017).

Regarding the comparison to tower albedo observations, it appears reasonable to suggest the metallic structure of the flux tower, located within the sensors field of view, may slightly increase the reflection of radiation towards the sensor compared to the canopy, particularly in the PAR. It is thus possible that these measurements slightly underestimate albedo. Beyond model structure, overestimation of albedo by both models could result from incorrect optical properties used for leaves, soil, or bark. There is some level of uncertainty in the parameters used across sites, even though we collected the best available data. These data can also vary significantly through time, as the humidity conditions are known to strongly influence the soil (and likely the bark) optical properties in the NIR.

5.3 Implications of including stems in microclimate models

Wood structures absorb smaller amounts of PAR than NIR (see Fig. 2), and it is expected that including stems in the radiative transfer simulation will have limited effect on fluxes of carbon and latent heat which are largely a function of PAR absorption. Radiation absorption by stems will more strongly influence sensible heat fluxes because the radiative forcing is a function of both PAR and NIR absorption. This influence on sensible heat and the storage of heat in stems has implications for the energy budget closure at the whole canopy level. The inclusion of stems in a microclimate model will also affect the amount of radiation reaching the soil, both in the NIR and PAR. Figure 2 shows that the absence of stems in the radiative transfer model nearly doubles the amount of radiation absorbed by the soil in the NIR and PAR compared to simulations including stems. This has significant implications for the energy balance at the soil-air interface. In CLM5.0, when the stems are considered as a canopy element absorbing radiation, the under-canopy wind speed is modified empirically to calculate the aerodynamic resistances to heat and moisture between the soil and the canopy air (Eq. 2.5.133 in Lawrence et al., 2018), such empirical adjustment is likely required following the decrease in radiative forcing on the soil caused by stem radiation absorption. We will discuss in further details the implications related to energy balance and energy budget in a subsequent paper.

6 Conclusions

We modified the Norman (1979) model to provide estimates of radiation absorption by wood that are directly computed from canopy structure information. Taking the 3D radiative transfer simulations as a benchmark, we showed that the modified Norman model accurately reproduces vertical profiles of absorbed radiation by sunlit leaves, shaded leaves, and wood in dense broadleaf forests. We showed that the modified Norman's model estimates are in closer agreement with albedo observations than the original model, and those estimates are in very close agreement with estimates from the 3D ray tracing radiative transfer model.

The model developments pertaining to radiation absorption by wood structures have implications for land surface models, whether they are using a big leaf or multilayer canopy scheme to compute mass and energy fluxes. Models such as the Community Land Model (CLM) combine leaf area index and stem area index into a plant area index to calculate radiative transfer and surface fluxes, but the radiative transfer does not distinguish the absorption by leaves and wood. When biomass heat storage is included in the CLM canopy energy balance, ad-hoc adjustments are used to partition the absorbed solar radiation to leaves and wood (Swenson et al., 2019). The modified Norman's model avoids relying on ad-hoc relations for calculating the radiation fluxes absorbed by wood structures, and these fluxes are resolved vertically. As shown here, absorption of radiation by wood is a non-negligible component of the surface energy budget, especially in the NIR. This will likely lead to improved calculations of heat storage in wood material, as well as improved heat fluxes and longwave radiation emission from stems and branches, which influence within canopy air temperatures. The level of agreement obtained here between the 1D and 3D model estimates suggests algorithmic implementation of the stem addition as an absorbing element is correct and can be used to benchmark modifications to other radiative transfer models that are more commonly used in land surface models, like the two-stream model of Dickinson (1983) and Sellers (1985). Such model developments, including further validation in needleleaf forests, mixed forests, and sparse canopies, may also play a significant role in improving the fit between land surface model simulations and flux tower observations by increasing model correctness and reducing reliance on model empirical parametrisation.