Simulating the 3D Distribution of Absorbed Shortwave Radiation in a Tree Crown: Comparison of Simplified and Monte Carlo Models

Urban trees are expected to alleviate severe thermal environments in summer. The absorbed shortwave radiation is a key variable to estimate the microclimatic effects of trees. In numerical simulation tools for the urban and built thermal environment, the difficulty lies in how to calculate the absorbed radiation in terms of accuracy and computational load. Here we investigated the difference between a radiative transfer model calculating scattering by the Monte Carlo method and a simpler model that accounts for scattering by modifying the optical depth. The calculations of absorbed radiation were performed under various conditions of leaf normal distribution and zenith angle of incident radiation, using tree models represented by voxels with different leaf area density values. For the photosynthetically active radiation, the relative difference between models was less than 10% for most cases. For the near‐infrared radiation, the difference was also small for a large canopy consisting of several trees. The difference was 30%–60% and 0%–20% for zenith angles of incident radiation of 0° and 60°, respectively, for small isolated trees with a crown height and width of less than 5 m, which was a significantly different condition from a uniform canopy. This large difference was related to the uniformity of the canopy and was explained by the transmittance from each voxel to the outside of the canopy. The present study enables us to evaluate the applicability of the simple model and correct the difference from the Monte Carlo model.


Introduction
Urban trees are expected to alleviate severe thermal environments in summer through solar shading and transpiration.Planting trees, however, does not always produce a cooling effect.Urban trees are expected to be under severe water conditions because of significant vapor pressure deficits and limited water supply and root space (Asawa et al., 2017;Gillner et al., 2017;Kagotani et al., 2013;Osone et al., 2014;Su et al., 2021); therefore, the Bowen ratio (ratio of sensible heat to latent heat) will tend to increase (Asawa & Fujiwara, 2020).Even under conditions with sufficient transpiration, sunlit leaves warm the air due to the dominant sensible heat flux (Cook et al., 1964;Knoerr & Gay, 1965), and thus sunlit vegetation surfaces should be kept away from (sufficiently higher than) the pedestrian level (Manickathan et al., 2018).The solar radiation absorbed by trees is a key variable to evaluate the effects of trees because such radiation is the source of sensible heat, latent heat, and longwave radiation from leaves.The leaf absorptance varies significantly according to the wavelength domain: ∼0.8 and ∼0.2 for photosynthetically active radiation (PAR) and near-infrared radiation (NIR), respectively.Despite the Abstract Urban trees are expected to alleviate severe thermal environments in summer.The absorbed shortwave radiation is a key variable to estimate the microclimatic effects of trees.In numerical simulation tools for the urban and built thermal environment, the difficulty lies in how to calculate the absorbed radiation in terms of accuracy and computational load.Here we investigated the difference between a radiative transfer model calculating scattering by the Monte Carlo method and a simpler model that accounts for scattering by modifying the optical depth.The calculations of absorbed radiation were performed under various conditions of leaf normal distribution and zenith angle of incident radiation, using tree models represented by voxels with different leaf area density values.For the photosynthetically active radiation, the relative difference between models was less than 10% for most cases.For the near-infrared radiation, the difference was also small for a large canopy consisting of several trees.The difference was 30%-60% and 0%-20% for zenith angles of incident radiation of 0° and 60°, respectively, for small isolated trees with a crown height and width of less than 5 m, which was a significantly different condition from a uniform canopy.This large difference was related to the uniformity of the canopy and was explained by the transmittance from each voxel to the outside of the canopy.The present study enables us to evaluate the applicability of the simple model and correct the difference from the Monte Carlo model.
Plain Language Summary Evaluation of the microclimatic effects of trees (solar shading and transpiration) is important for the designing of urban areas to alleviate severe thermal environments in summer.The absorbed solar radiation is a key variable to estimate the microclimatic effects of trees because the absorbed energy is divided into that used to heat the air and that used for transpiration.Numerical simulation is a feasible way to estimate the absorbed radiation.The Monte Carlo method is commonly used to simulate the absorbed radiation with consideration for the scattered radiation by leaves.However, current tools for simulating the thermal environment of urban and built spaces use a simpler method wherein the scattered radiation is implicitly considered by modifying a parameter (optical depth) in the radiative transfer equation.We here used three-dimensional tree data to compare the absorbed radiations simulated using the simple and Monte Carlo methods.For photosynthetically active radiation where leaf absorptance is high, the difference between models was small regardless of tree size.For near-infrared radiation where leaf absorptance is low, the difference was notable for small trees.We investigated how the difference between models was related to the tree's structural characteristics to clarify the applicability of the simple method.
• We compare the absorbed shortwave radiation between models with and without explicit treatments of leaf scattering using 3D tree models • The relative difference in absorbed near-infrared radiation between models exceeds 30% in many cases for small isolated trees • The large difference in individual voxels is due to the optical depth in surrounding directions, which represents the canopy uniformity

Supporting Information:
Supporting Information may be found in the online version of this article.
In the recent version of the ENVI-met, which is one of the most commonly used simulation models of urban thermal environment according to a review provided by Jänicke et al. (2021), a radiation transfer module that accounts for the scattered radiation by leaves in an approximate manner was incorporated (Simon et al., 2020).
The module is based on Goudriaan's simple model (Goudriaan, 1977), in which the decrease of net radiation at a certain wavelength domain through a uniform canopy is represented by modifications to the optical depth using the Beer-Lambert law.This simple module can be easily incorporated into simulation models of the urban thermal environment because the calculation method is similar to that of radiation extinction through a canopy.
This study sought to elucidate the applicability of Goudriaan's simple model to the calculation of solar radiation absorbed by trees with a nonuniform crown structure.To represent the crown structure, light detection and ranging (LiDAR)-derived voxel data of leaf area density (LAD: total one-sided leaf area per unit volume) were used.
The absorbed PAR (APAR) and ANIR calculated using a voxel-based simplified radiative transfer model based on Goudriaan's simple model (VG model) were compared with those calculated by an MC model.To elucidate the types of tree canopy to which the VG model can be applied, we investigated how the difference in absorbed radiation between models is related to canopy structural characteristics.The details of the radiative transfer models are described in Section 2. The LAD data and simulation conditions are provided in Section 3. The difference between models and the structural characteristics explaining the difference are reported in Section 4.

VG Model
The solar radiation (PAR or NIR) absorbed by a voxel was calculated on the basis of Goudriaan's simple model as:  (Kobayashi, 2012), in which the absorbed radiation per voxel is calculated by a voxel center approach (i.e., multiplying the absorbed radiation per leaf area at the center of the voxel by LAD).The optical depth is multiplied by  √  based on Goudriaan's simple model, by which the decrease of net radiation (useable radiation within the canopy, including radiation passing through gaps and radiation scattered by leaves) is represented.This representation is theoretically correct for a canopy where horizontal leaves with equal reflectance and transmittance are randomly distributed (Cowan, 1968) and is a good approximation for inclined leaves (Goudriaan, 1977).We added two factors to the original equation of FLiESvox: the canopy reflectance ρ c (μ,θ b ) and the correction factor for the voxel center approach f err .
Canopy reflectance is required to calculate the net radiation.We adopted the following widely used equation (Spitters, 1986) in the present study: where a 1 and a 2 are coefficients.This equation describes a directional-hemispherical reflectance of a uniform canopy.For spherical leaf normal distribution (LND), values of 2 and 1.6 have been used for a 1 and a 2 , respectively (Spitters, 1986).To apply the model to other types of LND, coefficients for erectophile, plagiophile, and planophile LNDs (Bunnik, 1978) were determined using the following MC model.Details are described in Appendix A. The coefficient values are shown in Table 1.The spherical LND was included in the calculation, and the resulting values differed from those obtained previously.The values shown in Table 1 are for a canopy with a ground reflectance of zero.The reason for this treatment is that the reflected radiation can be calculated separately.
A correction factor was introduced to account for errors due to the voxel center approach.The voxel center approach is a good approximation when a balance exists between the decrease in radiation reaching the voxel center with increasing LAD and the increase in absorptivity of the voxel with increasing LAD.However, the balance shifts with increasing LAD, leading to an underestimation of absorbed radiation.To adjust for the underestimation, we obtained a correction factor by comparing the absorbed radiation calculated using the voxel center approach and that calculated considering all possible paths through a voxel.Details are described in Appendix B. The obtained equation is where   ′ vox is the modified optical depth of the target voxel (i = 1 in Equation 2).

MC Model
Our MC radiative transfer model was based on the FLiES (Kobayashi & Iwabuchi, 2008) top of the canopy was set arbitrarily).We confirmed the validity of our simulation code by comparing its results with those obtained using the FLiES model (Text S1, Figures S1 and S2, Tables S1 and S2 in Supporting Information S1).Note that a voxel size of 0.05 m was used in this comparison to reduce the difference in tree shape between models (geometric shape in FLiES and voxel representation in our MC model).After confirming that the photon paths were properly traced using our MC model by means of the above-described comparison, our MC model was applied to the LiDAR-derived tree models with larger voxel sizes shown in Section 3.1.

LAD Data
Voxel-based LAD data sets for five different isolated trees were used: Cerasus × yedoensis, Cinnamomum camphora, Ginkgo biloba, Quercus myrsinifolia, and Zelkova serrata.These trees were planted in an experimental field with an open area of 8,800 m 2 in the city of Miyoshi, Aichi Prefecture, Japan (35.1355°N, 137.1001°E).
We selected these species in part because they have different crown shapes and structures.Each tree was measured by a terrestrial LiDAR system from four positions surrounding the tree.The LiDAR measurements were conducted on 29 August 2012, and the target trees were 13 years old at that time.The LAD with a voxel size of 0.3 m was estimated on the basis of contact frequency between laser beams and leaves.Details of the LiDAR measurements and LAD estimation are described in Oshio and Asawa (2020).LAD represents the volumetric density of leaves; for example, a forest canopy with an LAI of 5 m 2 m −2 and a height of 10 m has an LAD of 0.5 m 2 m −3 .In the present study, LAD was given for each voxel to represent the 3D structure of the tree crown (i.e., the LAD of each voxel was multiplied by the voxel height and integrated in the vertical direction, yielding LAI). Figure 1a shows photographs and LAD data of the target trees.The tree height and LAI are also shown in the figure.The LAI was obtained using LAD data as the ratio of the total one-sided leaf area within the crown to the projected area of the crown (projected area of voxels with LAD > 0 on the horizontal surface).
As mentioned in previous sections, the modification of the optical depth as in Equation 2and canopy reflectance as in Equation 3 are based on the radiative transfer theory for a uniform canopy; therefore, the applicability of these equations is expected to depend on the canopy size.Then, for comparison, trees with larger crowns were also used.Airborne LiDAR-derived LAD data (Oshio & Asawa, 2016;Oshio et al., 2015) was employed.The data covers a part of Hisaya-Odori Street in the city of Nagoya, Aichi Prefecture, Japan (35.1709° N, 136.9085°E).
The LAD was estimated using multi-return airborne LiDAR data with a density of 20 points m −2 obtained on 6 September 2010.The voxel size was 0.5 m.Details of the airborne LiDAR observation and LAD estimation are described in Oshio et al. (2015).From this data, we selected two trees: an isolated Z. serrata tree and a C. camphora canopy.Figure 1b shows photographs and LAD data of these trees.The C. camphora canopy was a row of trees in the north-south direction.
Woody elements were not considered in the simulation.In short, the present simulation was focused on investigating characteristics of radiative transfer in a canopy with some degree of dense foliage.

Simulation Conditions
Simulations were conducted under various conditions of the LND, zenith angle of incident radiation beam (BZA), and wavelength domain (Table 2).For the LND, four typical distributions were considered: erectophile, spherical, plagiophile, and planophile.Note that the species-specific LND was used to calculate the LAD (Oshio & Asawa, 2020;Oshio et al., 2015).Here, these four distributions were used to investigate the influence of LND on the difference between the VG and MC models, irrespective of the actual distribution.For BZA, three different settings were used with a constant azimuth angle of 0° (east direction).To calculate the absorbed radiation derived from direct and sky-diffuse solar radiations using the VG model, a ray is traced in the direction of the sun, and rays are traced in other directions, respectively.Each of the above settings of incident angles could be regarded as corresponding to one of them (i.e., the solar direction or one other direction in the sky).Whether the direct or sky-diffuse component was not specified because we just wanted to investigate the difference between models depending on the zenith angle.The above setting of azimuth angle had no special meaning since there were no special directional characteristics for isolated trees.For the C. camphora canopy, the radiation beam entered from the perpendicular direction of the row.For the wavelength domain, PAR and NIR were considered.The leaf absorptance was set to 0.8 and 0.2 for PAR and NIR, respectively.The leaf reflectance and transmittance were set to half of one minus absorptance.These settings are often used as typical values for leaves (Kobayashi et al., 2002;Ma et al., 2018) and were also close to the actual values of the target trees in the present study (Text S2, Figure S3, Tables S3 and S4 in Supporting Information S1).Since the objective of this study was not to evaluate the absorbed radiation of particular trees but to examine the difference between radiative transfer models in general, typical values were used.The radiation flux on a plane perpendicular to the beam direction was set to 400 W m −2 for both domains.The reflectance of the ground surface was set to zero to evaluate the absorbed 10.1029/2023JD038612 6 of 22 radiation purely derived from the downward radiation incident on the canopy (radiation reflected at the ground can be calculated separately).The isolated boundary condition was used (tracing was terminated for rays and photons escaping the canopy).For MC simulation, the number of photons was 25,000 per unit voxel area (the square of the length of one side of the voxel).For the above conditions, APAR and ANIR were calculated using the VG model (APAR VG and ANIR VG ) and MC model (APAR MC and ANIR MC ).

Evaluation of the Difference Between Models
To evaluate the difference in APAR and ANIR between the VG and MC models, the relative difference was calculated: where ∆ r is the relative difference for PAR or NIR [%],  tot_VG is the total APAR VG or ANIR VG in the tree [W], and  tot_MC is the total APAR MC or ANIR MC in the tree [W].To evaluate the variation in the difference between VG and MC models, the normalized standard deviation (NSD) was calculated: where s nor is the NSD for PAR or NIR [%], s is the standard deviation of the difference between APAR VG and APAR MC or that between ANIR VG and ANIR MC in the tree [W m −3 ], and  MC is the mean value of APAR MC or ANIR MC in the tree [W m −3 ].In addition, to evaluate the correlation between the VG and MC models, the coefficient of determination of linear regression (R 2 ) was calculated for PAR using APAR VG and APAR MC and for NIR using ANIR VG and ANIR MC .
The amount of absorbed radiation of a voxel varies with LAD, the G function value (i.e., LND and BZA), and the incident radiation.In order to evaluate how the difference in absorbed radiation between models is related to canopy structural characteristics without being affected by such conditions, the normalized difference was calculated as where ∆ nor is the normalized difference

APAR
The vertical section of APAR is shown in  LAD; therefore, the small difference between APAR VG and APAR MC for voxels consisting of a crown surface on the side where solar radiation entered (upper and east surfaces for the BZA of 0° and 60°, respectively) is natural.Good agreement was also obtained in the inner part of the crown.It was expected that the difference between APAR VG and APAR MC was smaller for C. camphora on Hisaya-Odori Street because its large canopy was closer to the condition of a uniform canopy; however, in the figure, there is little apparent difference even for the small isolated Q. myrsinifolia tree.
Figure 3 shows the scatter plot of APAR VG and APAR MC .In addition to the results for planophile LND, those for erectophile LND, which exhibited the smallest relative difference within the results for NIR, are shown.The scatter plots clearly show the agreement between models.Although APAR VG tended to be greater than APAR MC for voxels having large APAR value, especially for planophile LND, most data were close to the 1:1 line.
The relative difference, NSD, and R 2 for erectophile and planophile LNDs are shown in Figure 4 (those for spherical and plagiophile LNDs are shown in Figure S4 in Supporting Information S1).Although the relative difference for isolated trees was greater than that for C. camphora on Hisaya-Odori Street, it was less than 10% for most cases.The relative difference varied with BZA, and the variation was greater for planophile LND than for erectophile LND.Even for PAR, where the leaf absorptance is high, LND and BZA affected the results.
The G function values for BZAs of 0°, 30°, and 60° were 0.42, 0.45, and 0.5, respectively, for erectophile LND and 0.84, 0.72, and 0.47, respectively, for planophile LND.For trees in the experimental field, although the G function value was small for erectophile LND with a BZA of 0°, the probability of the radiation beam contacting a leaf was not very low due to the vertically long crown shape.The G function value increases slightly for erectophile LND with BZA 60°; however, the contact probability did not increase much due to the short path length through the crown.For planophile LND, the contact probability was high when BZA was 0° due to the large G function value and long path length.In this case, the amount of the scattered radiation escaping from the crown was expected to be large because there were a smaller number of surrounding voxels containing leaves due to the vertically long crown shape.Therefore, the calculation based on Goudriaan's simple model assuming a uniform canopy overestimated the APAR, especially for planophile LND with a BZA of 0° (Figure 4d).
The NSD was less than 10% in many cases for erectophile LND; however, greater NSD and variation in NSD with BZA were seen for planophile LND (Figures 4b and 4e).The correlation was very high (Figures 4c and 4f), although a certain level of high correlation was natural because the absorbed radiation depended on the LAD of voxels.The greater NSD and smaller R 2 for G. biloba could be attributed to the heterogeneous foliage distribution.In summary, the difference between APAR VG and APAR MC was small even for small isolated trees: the relative difference and NSD were less than 10% and 20%, respectively, for most cases.

ANIR
The vertical section of ANIR is shown in Figure 5. Unlike APAR, ANIR was relatively large in the inner part of the crown.Although the distribution of ANIR VG was similar to that of ANIR MC , the difference in the absolute value can be read from the figure, especially for Q. myrsinifolia in the experimental field with a BZA of 0°.
Figure 6 shows the scatter plot of ANIR VG and ANIR MC for erectophile and planophile LNDs.Overall, the data variability is greater for NIR than for PAR.For Q. myrsinifolia with a BZA of 0°, data are clearly plotted above the 1:1 line regardless of LND.For C. camphora on Hisaya-Odori Street, most data were close to the 1:1 line; however, ANIR VG was clearly greater than ANIR MC for voxels with large ANIR, for planophile LND with a BZA of 0°.
The relative difference, NSD, and R 2 for erectophile and planophile LNDs are shown in Figure 7 (those for spherical and plagiophile LNDs are shown in Figure S5 in Supporting Information S1).The relative difference   was significantly greater than that for APAR (the range of the vertical axis is different from Figure 4).For small trees in the experimental field, the relative difference varied with BZA, and the variation for planophile LND was greater than that for erectophile LND.Although these characteristics can be explained in the same way as PAR, the difference between ANIR VG and ANIR MC was quite large because a large amount of scattered radiation appeared to escape from the crown due to the low leaf absorptance and the small crown size.The VG model overestimated the scattered radiation that can be used in the crown.For erectophile LND, the dependence of the relative difference on BZA was small for PAR but increased for NIR (Figures 4a and 7a).The influence of crown shape seemed to become more pronounced: the upper part of the crown on which the radiation beam with a BZA of 0° entered was far from the assumption of a uniform canopy due to the vertically long shape.For C. camphora on Hisaya-Odori Street, which had a large canopy, the relative difference was much smaller than that for small trees and comparable with that for the APAR of small trees.The NSD was also significantly greater for ANIR than for APAR, especially for small trees (Figures 7b and 7e).The NSD showed characteristics similar to those for the relative difference.The R 2 for ANIR was lower than that for APAR but still high (around 0.96) (Figures 7c and 7f).
In summary, C. camphora on Hisaya-Odori Street appeared to act as a uniform canopy, and the ANIR was well estimated under various conditions by the VG model.In contrast, for small trees, the difference between ANIR VG and ANIR MC was very large in some cases and depended on BZA: the relative difference was 40%-60% with an NSD of 50%-80% for 0° BZA and 0%-10% with NSD of 20% for 60° BZA, under conditions with planophile LND.

Analysis of the Difference in ANIR Between Models
As shown in the previous section, the difference in ANIR between models varied among trees and showed even greater variation by canopy size.Therefore, we investigated which variables represent canopy structural characteristics that affect differences between the models, in order to quantitatively clarify the applicability of the VG method.

Variation With Height
The variations in the normalized difference (Equation 7) with height for erectophile and planophile LNDs are shown in Figures 8 and 9, respectively (those for spherical and plagiophile LNDs are shown in Figures S6 and S7 in Supporting Information S1, respectively).The absolute value of the normalized difference can be extremely large when LAD is close to zero (i.e., division by a value close to zero, see Equation 7); therefore, data with LAD > 0.002 m 2 m −3 were used in Figures 8 and 9.In Figure 8, for small trees in the experimental field, the normalized difference decreased toward the lower part of the crown.The difference between BZAs was small.These results indicate that the normalized difference is related to the optical depth.More specifically, when BZA was 0°, a long geometric distance is required for the incident radiation beam to reach the lower part of the crown due to the vertically long shape of the crown, but the small G function value yielded a moderate optical depth.
When BZA was 60°, the geometric distance was short, but the G function value was large, yielding an optical depth comparable to that for BZA 0°.The effects of BZA, LND, and crown shape appeared to be balanced, and similar characteristics were obtained for different BZAs.Although the difference in results between small trees was not remarkable, a somewhat greater normalized difference at the upper part of the crown (interception of the linear regression line) for G. biloba and Q. myrsinifolia can be attributed to the pointed tree tops, which form a canopy that is far from uniform.For trees on Hisaya-Odori Street, the variation in the normalized difference with height was small.The large isolated Z. serrata tree showed the same characteristics as the C. camphora canopy.At the top of the crown (topmost voxel), the normalized difference for trees on Hisaya-Odori Street was comparable with that for small trees (∼0.2).However, the normalized difference for trees on Hisaya-Odori Street promptly decreased to a value similar to that at the lower part of the crown.It appeared that not only the optical depth in the direction of the incident radiation beam but also the amount of surrounding foliage affected the difference between models.
Regarding planophile LND (Figure 9), for small trees in the experimental field, the difference between BZAs was notable, which was different from the results in Figure 8. Greater interception and negative slope were obtained for 0° BZA.As described in Section 4.1, the G function value was large when BZA was 0°, and the probability of an incident radiation beam contacting the leaves and scattering was high.The VG model did not account for scattered radiation that escaped from the canopy and became unavailable.Meanwhile, the optical depth increased toward the lower part of the crown, and approached characteristics of a uniform canopy, yielding a small normalized difference.This effect was observed even for trees on Hisaya-Odori Street.

Relationship Between Model Difference and Canopy Structural Variables
According to the previous section, the optical depth appeared to be related to the normalized difference.Therefore, the transmittance in the direction of the incident radiation beam (the exponential part of Equation 1) was used: The presence of surrounding foliage also appeared to be important.The mean of transmittance in all directions was calculated as where T mean is the mean transmittance [−], N z and N a are the numbers of zenith and azimuth angles for which transmittance is calculated [−], τ′(θ i ,φ j ) is the modified optical depth to the direction of ith zenith angle θ i and jth azimuth angle φ j [−].The calculation was conducted with an interval of 5° (N z = 36 and N a = 72).
In addition, weighted transmittance, which is a variable considering characteristics of the above two variables was calculated as  ).We also calculated  _down that is the weighted transmittance in the downward direction (the same calculation as for  _up was performed for negative zenith angles).The calculations of  _up and  _down were conducted separately to take into account the difference between the upper and lower parts of the crown (at the upper part of the crown,  _up was expected to be high and  _down was expected to be low; at the lower part of the crown, the opposite was true).The relationship between the normalized difference and variables is shown in Figure 10.As in the analyses in Figures 8  and 9, data with LAD > 0.002 m 2 m −3 were used.In the figure, data are plotted without distinguishing between BZAs because the difference between BZAs was expected to be explained by variables.The normalized difference increased with increasing transmittance in the beam direction (Figures 10a-10d).Q. myrsinifolia showed a greater normalized difference than C. camphora on Hisaya-Odori Street, even at the same transmittance.The variation in the normalized difference at a transmittance value was larger for C. camphora on Hisaya-Odori Street than for Q. myrsinifolia.
When the mean transmittance was considered (Figures 10e-10h), the normalized difference increased with increasing mean transmittance when the transmittance in the beam direction was high and hardly depended on mean transmittance when the transmittance in the beam direction was low.The histogram in the figure shows that the fraction of voxels having high transmittance in the beam direction was high for Q. myrsinifolia.Meanwhile, there were voxels with low transmittance in the beam direction due to the small but dense crown.For the mean transmittance, the difference between trees was clear: few voxels had mean transmittance lower than 0.5 for Q. myrsinifolia, but the peak of the distribution was about 0.5 for C. camphora on Hisaya-Odori Street.
Figures 10i-10l shows that the normalized difference increased with a simultaneous increase in both weighted transmittances.Considering the result that the normalized difference was large at the upper part of the crown (Figures 8 and 9), one might have expected that the normalized difference would increase with an increase in the weighted transmittance in the upward direction and a decrease in that in the downward direction.However, this was not the case, because the more leaves in the downward direction, the closer the canopy characteristics approach those of a uniform canopy.These characteristics could not be explained by simple variables, such as height and transmittance in the beam direction.The range of transmittances differed between trees, and few voxels of Q. myrsinifolia had low transmittances.

Empirical Correction for Model Difference
The equation estimating the normalized difference was obtained from the relationship shown in Figure 10 in order to investigate the generality of the relationship and enable corrections when the difference between models was large.Three equations were obtained using linear regression: Since ANIR is a positive quantity, the following treatment was used when calculating the sum of squares of the residuals during the regression.If the corrected ANIR was negative using the estimated normalized difference value, the estimate of the normalized difference was changed to a value such that the corrected ANIR was zero, and the residual was calculated using the changed value.The sum of squares of the residuals was calculated while changing the parameter values by 0.001, and parameters offering the minimum sum value were used as the optimal coefficients.Coefficients were obtained for each tree and LND (Table S5 in Supporting Information S1).In the following analysis, coefficients obtained from Cerasus × yedoensis with spherical LND (Table 3), which were the intermediate values of the small trees in the experimental field, were applied to other trees and LNDs to test the applicability of the coefficients.During the correction, if the corrected ANIR was negative, it was changed to zero, as in the regression calculations.
The scatter plot of the corrected ANIR VG and ANIR MC is shown in Figure 11.The results shown are for Q. myrsinifolia in the experimental field and C. camphora on Hisaya-Odori Street and for erectophile and planophile LNDs, which are different trees and LNDs from those used to obtain the regression coefficients.When the correction was conducted using the transmittance in the beam direction, data were plotted closer to the 1:1 line than before the correction for Q. myrsinifolia in the experimental field.In contrast, most data were below the 1:1 line after correction for C. camphora on Hisaya-Odori Street.Some voxels showed ANIR VG of zero but ANIR MC greater than zero due to the treatment of negative ANIR in the correction; however, the number of such voxels was small.When the correction was conducted using the weighted transmittances in the upward and downward directions, the results changed slightly for Q. myrsinifolia in the experimental field, although the variation became larger for 0° BZA, and data were closer to the 1:1 line for 60° BZA.For C. camphora on Hisaya-Odori Street, the results were improved, except for 0° BZA with planophile LND.The difference in results between before and after correction was small (Figures 6e-6h and 11m-11p).
The relative difference, NSD, and R 2 for the corrected ANIR VG for erectophile and planophile LNDs are shown in Figure 12 (results for spherical and plagiophile LNDs are shown in Figure S8 in Supporting Information S1).
When the correction was conducted using the transmittance in the beam direction, the relative difference was less than 10% in most cases for small trees in the experimental field, although there was still some dependence of the difference on BZA.The same was true for Z. serrata on Hisaya-Odori Street.However, the relative difference was −20% to −30% for C. camphora on Hisaya-Odori Street.These results indicate that the transmittance in the beam direction can be used for corrections for isolated trees whose crowns are not very large, but not for explaining the relative difference in general.When the correction was conducted using the transmittance in the beam direction and mean transmittance, the relative difference for C. camphora on Hisaya-Odori Street was reduced.However, the relative difference was large for small trees with planophile LND.This indicates that the regression coefficients obtained from Cerasus × yedoensis with spherical LND just happened to be appropriate for C. camphora on Hisaya-Odori Street with planophile LND, but the transmittance in the beam direction and mean transmittance did not fully account for the crown characteristics.When using these explanatory variables, the multiple correlation coefficient was high in the regression of each tree (Table S5 in Supporting Information S1).However, it was difficult to obtain general coefficients.When the correction was conducted using the weighted transmittances, the relative difference was less than 10% in most cases, including for C. camphora on Hisaya-Odori Street.The weighted transmittances appeared to represent the crown structural characteristics that were related to the difference in ANIR between models.Comparing Figures 7b, 7e, 12n, and 12q, it can be seen that the NSD was not decreased by the correction.However, the relative difference close to zero and the small difference between trees is important, because it is not often assumed that one would like to know the ANIR at a specific voxel.

Table 3
Coefficients for the Estimation of Normalized Difference

Discussion
The simple radiative transfer model developed by Goudriaan (1977) has been compared to models explicitly considering scattering.Indeed, Goudriaan (1977) conducted such a comparison upon introducing the simple   2020) also compared the radiation flux calculated using the tree model with LAD distribution with that calculated using the tree model with a uniform LAD (the same LAD for all voxels).As is already well known-and in fact was the motivation for our use of the voxel-based tree model with LAD distribution-the tree model with a uniform LAD showed significantly smaller radiation flux in the tree crown.Tree models with LAD distribution are expected to become more important for simulating the thermal environment of urban and built spaces.This is because such models can bridge the gap between the design of urban and built spaces utilizing trees and the impact of trees on the thermal environment.Using computer-generated tree models is useful in predicting the thermal environment, but generating tree models from LiDAR measurements is also useful for the assessment of physical spaces.LAD data with voxel sizes of 0.3 and 0.5 m were used in the present study since such data can be obtained using LiDAR measurements and used in numerical simulation.In fact, the difference in the calculated absorbed radiation between VG and MC models was not highly dependent on the voxel size (Text S3 and Figure S9 in Supporting Information S1).Rather, the larger the voxel size, the larger was the calculated value of absorbed radiation itself.Naturally, as the voxel size increased, the vertical distribution of absorbed radiation could not be represented (Figures S10 and S11 in Supporting Information S1).If the voxel size was small, even the VG model could represent the shape of the vertical distribution.Therefore, when considering where solar radiation is absorbed and where sensible heat is emitted, the use of LAD data with a voxel size similar to that used here is recommended.The present study shows that ANIR obtained using the simple model is significantly greater than that obtained using the MC model even when a tree model with a small voxel size is used.Incidentally, as for correction, coefficients obtained from small voxel sizes could be used for larger voxel sizes since variables depending on voxel size were not used (Figure S9 in Supporting Information S1).In general, data of 3D LAD distribution have not been readily available.However, airborne high-density point clouds are now becoming publicly available (e.g., Shizuoka point cloud DB, 2022; USGS 3D elevation program, 2022), and methods for estimating the 3D distribution of LAD or plant area density using airborne point clouds have been developed (Arnqvist et al., 2020;de Almeida et al., 2019;Halubok et al., 2022;Oshio et al., 2015).In addition, small and inexpensive laser scanners have been developed, and point clouds are readily available by measurements from the ground or an unmanned aerial vehicle (Hu et al., 2021;Lu et al., 2020).
Our results indicate that the simple method presented herein with 3D LAD data is a useful approach to more accurately account for scattering within the framework of radiation transfer calculations in conventional simulation tools for urban and built thermal environments.For the correction described in Section 4.3.3,coefficients obtained from one species could be applied to trees with different species and sizes, demonstrating the applicability of coefficients.However, the applicability of coefficients to other tree types such as coniferous trees needs further investigation.Of course, these corrections are not necessary if an MC model is used, but the present study can help in selecting which of the various models with different accuracies and computational loads to employ.For computational load, the calculation speed of the MC model depends on the number of photons.For the trees used in the present study, the MC model took longer computation time than the VG model when calculated with the number of photons required to obtain stable results (Text S4, Figure S12, and Table S6 in Supporting Information S1).Comparison in a larger study area including buildings is needed to fully investigate the difference in computational load between models.
Radiative transfer models that can perform a calculation using tree models with individual leaves and branches or those with very small voxels (∼0.1 m) have been developed (Li et al., 2018;Widlowski et al., 2014).
Our tree models were represented by voxels with a spatial scale of foliage.This representation is a realistic treatment considering that in addition to the computational load, the evaluation of the thermal environment involves computational fluid dynamics simulations as well as radiative transfer calculations (Grylls & van Reeuwijk, 2021;Manickathan et al., 2018;Oshio et al., 2021).Woody elements were not considered in the present study.Talbot and Dupraz (2012) indicated that the influence of woody area density on the PAR transmittance under a leafy tree canopy was small.According to Kobayashi et al. (2012), the radiation absorbed by woody elements accounts for a notable proportion of the radiation absorbed by the savanna canopy where the LAI is low, especially for NIR.Widlowski et al. (2014) also showed that omitting elements resulted in a non-negligible bias in the calculation of the bidirectional reflectance factor for the savanna canopy.For denser canopies, Ebengo et al. (2021) compared the hyperspectral image obtained by airborne observation and that obtained by numerical simulation in rainforests and showed that the difference in bias between treatments of woody elements was small even for NIR.Schneider et al. (2014) also showed that the simulated at-sensor radiance in the NIR domain decreased only slightly by adding woody elements in the calculation for an old temperate mixed forest.Therefore, the impact of neglecting woody elements on the absorbed radiation is expected to be small for leafy trees such as those used in the present study.This was partially confirmed quantitatively for the tree in the experimental field (Text S5 in Supporting Information S1).In the NIR region, the influence of woody elements was greater than in the PAR region, but the ANIR by woody elements was expected to be 10%-15% of the total (Table S7 in Supporting Information S1), not enough to change the vertical distribution of the absorbed radiation (Figures S13 and S14 in Supporting Information S1).Practically, however, it is important to take into account the decrease in leaf density due to seasonal changes, tree health, and pruning.For the tree in the experimental field, the ANIR by woody elements exceeded 30% of the total when LAI dropped to 1 (Figure S15 in Supporting Information S1).In the future, methods of representing trees including woody elements and calculating radiation transfer will be investigated with an eye toward a balance between accuracy and practicality.

Conclusions
We used urban tree models represented by voxels having different LAD values to compare the absorbed radiation calculated using a simple method based on Goudriaan's simple model and that calculated using an MC method.The difference in APAR between models was less than 10% for most conditions, even for small isolated trees.The difference in ANIR between models was also small for the large C. camphora canopy, which showed similar characteristics to a uniform canopy.For small isolated trees, the difference in ANIR between models was large and varied with the BZA: 30%-60% and 0%-20% for the BZAs of 0° and 60°, respectively.The large difference in ANIR was related to the uniformity of the canopy and was explained by the transmittance from each voxel to the outside of the canopy.In particular, assuming a plane facing the direction of the beam, the transmission in each direction weighted by the view factor from that plane corresponded well with the difference between models.According to the present results, the equation relating the difference in ANIR between models to the weighted transmittances, which was obtained from a small Cerasus × yedoensis tree with spherical LND, can be applied to various trees with different sizes and LNDs.
Our results provide a quantitative guide to the validity of the method using voxel-based tree models and Goudriaan's simple model, which is implemented in current simulation tools for the urban and built thermal environment.It is hoped that the present study will lead to a correct understanding of the impact of urban trees on the thermal environment.
R abs is absorbed radiation per voxel [W m −3 ], I b is the radiation flux on a plane perpendicular to the radiation beam at the top of the canopy [W m −2 ], ρ c (μ,θ b ) is the canopy directional-hemispherical reflectance [−], μ is the leaf absorptance [−] (the absorptance is the radiant energy absorbed by the target object divided by the incident radiant energy, i.e., unitless, and this is also true for reflectance and transmittance), θ b and φ b are the zenith and azimuth angles of the incident radiation beam at the top of the canopy, respectively [°], G(θ b ) is the mean projection area of a unit leaf area onto a plane perpendicular to the beam direction [−], u is the LAD of the target voxel [m 2 m −3 ], f err is the correction factor for the voxel center approach [−], τ′(θ b ,φ b ) is the modified optical depth (the optical depth is multiplied by  √  ) to the direction of the incident radiation beam [−], N is the number of voxels through which the radiation beam passes before reaching the center of the target voxel (i = 1 for the target voxel) [−], u i is the LAD of the i th voxel, and l i is the path length of the i th voxel [m].Equation 1 is based on a previously developed voxel-based model of a forest light environment simulator, FLiESvox

Figure 1 .
Figure 1.Photographs and voxel models of trees used in the present study: (a) trees in the experimental field; (b) trees on Hisaya-Odori street.The tree height (H) and LAI are shown in each photograph.
Figure 2. Results for Q. myrsinifolia with planophile LND are shown in the top row.As shown in the next section, the relative difference was greater for NIR than for PAR.Within the results for NIR, Q. myrsinifolia with planophile LND showed the greatest relative difference; therefore, these results are shown in the figures in the next section.For PAR, results for the same species and LND are shown for consistency.For purposes of comparison, results for C. camphora on Hisaya-Odori Street are shown on the bottom row.As shown in Figure 2, the absolute value and spatial distribution of APAR VG agreed very well with those of APAR MC for both trees.If the incident radiation on a voxel is the same, the absorbed radiation varies with Parameter Value or setting LND Erectophile, Spherical, Plagiophile, Planophile BZA 0°, 30°, 60°W avelength domain (leaf absorptance) PAR (0.8), NIR (0.2)

Figure 2 .
Figure 2. Vertical section of LAD, APAR MC APAR VG , and the difference: (a) Q. myrsinifolia in the experimental field; (b) C. camphora on Hisaya-Odori Street.For APAR, results with planophile leaf normal distribution for BZAs of 0° and 60° are shown.

Figure 8 .
Figure 8. Variation in the normalized difference (Δ nor ) for ANIR with erectophile leaf normal distribution according to the distance from the tree top for (a-e) trees in the experimental field and (f, g) trees on Hisaya-Odori Street.The box-and-whisker plot shows Δ nor for each height level (minimum, lower quartile, median, upper quartile, and maximum).The mean value (filled circle) and linear regression line to the mean value (solid line) are also plotted.Results for BZAs of 0° and 60° (blue and red, respectively) are plotted side by side with slight shifts (e.g., in the case of trees in the experimental field, the results for the top voxel are plotted at 0.075 and 0.225 m for BZAs of 0° and 60°, respectively).

Figure 9 .
Figure 9. Similar to Figure 8 but for planophile leaf normal distribution.

Figure 10 .
Figure10.Variation in the normalized difference (Δ nor ) for ANIR with (a-d) the transmittance in the beam direction (T beam ), (e-h) T beam and the mean transmittance (T mean ), and (i-l) the weighted transmittances in the upward and downward directions (T w_up and T w_down , respectively).Results for different trees and LNDs are shown: (a, b, e, f, i, j) Q. myrsinifolia in the experimental field;(c, d, g, h, k, l)  C. camphora on Hisaya-Odori Street; (a, c, e, g, i, k) erectophile LND; (b, d, f, h, j, l) planophile leaf normal distribution.The bar chart displayed along the axis shows the fraction of number of data in each bin (the vertical axis ranges from 0 to 0.1).

Figure 11 .
Figure 11.Scatter plot of ANIR MC and the corrected ANIR VG (ANIR VG_cor ): (a-h) correction using the transmittance in the beam direction (T beam ); (i-p) correction using the weighted transmittances in the upward and downward directions (T w_up and T w_down , respectively).Results for different trees and LNDs are shown: (a-d, i-l) Q. myrsinifolia in the experimental field; (e-h, m-p) C. camphora on Hisaya-Odori Street; (a, b, e, f, i, j, m, n) erectophile LND; (c, d, g, h, k, l, o, p) planophile leaf normal distribution.

Figure 12 .
Figure 12.Variation in(a, d, g, j, m, p)  relative difference (RD), (b, e, h, k, n, q) normalized standard deviation, and (c, f, i, l, o, r) R 2 for ANIR with BZA: (a-f) correction using the transmittance in the beam direction (T beam ); (g-l) correction using T beam and the mean transmittance (T mean ); (m-r) correction using the weighted transmittances in the upward and downward directions (T w_up and T w_down , respectively).Results for (a-c, g-i, m-o) erectophile leaf normal distribution (LND) and (d-f, j-l, p-r) planophile LND are shown.
. The individual trees are represented by geometric shapes (cone, cylinder, etc.), not voxels, in the FLiES model; therefore, a program code executing the same photon tracing as FLiES in voxel space was made.The atmospheric radiative transfer module is included in the FLiES model but was not included in the present study (the incident radiation on the

Table 2
Parameters Varied in the Examination are the ith zenith and jth azimuth angles in the coordinate system with the incident beam direction as the zenith [°].The calculation was conducted with an interval of 5° only for positive zenith angles ( where  _up is the weighted transmittance in the upward direction [−], S w is the sum of weights [−],