Advances in Modelling Radiative Transfer, Heat Storage and Turbulent Transport to Evaluate CO(2), Heat and Water Fluxes Over Broad-Leaved Forests: The CanVeg2 Model.

Advances in Modelling Radiative Transfer, Heat Storage and Turbulent Transport to Evaluate CO2 , Heat and Water Fluxes Over Broad‐Leaved Forests: The CanVeg2 Model - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Glob Chang Biol . 2026 Apr 18;32(4):e70867. doi: 10.1111/gcb.70867 Search in PMC Search in PubMed View in NLM Catalog Add to search Advances in Modelling Radiative Transfer, Heat Storage and Turbulent Transport to Evaluate CO 2 , Heat and Water Fluxes Over Broad‐Leaved Forests: The CanVeg2 Model Martin Béland Martin Béland 1 Digital Forest Lab, Department of Geomatics Sciences, Laval University, Quebec City, Quebec, Canada Find articles by Martin Béland 1, ✉ , Gordon B Bonan Gordon B Bonan 2 NSF National Center for Atmospheric Research, Boulder, Colorado, USA Find articles by Gordon B Bonan 2 , Tilden P Meyers Tilden P Meyers 3 NOAA/Air Resources Laboratory, Boulder, Colorado, USA Find articles by Tilden P Meyers 3 , J William Munger J William Munger 4 School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts, USA Find articles by J William Munger 4 , Hideki Kobayashi Hideki Kobayashi 5 Institute of Arctic Climate and Environment Research, Japan Agency for Marine‐Earth Science and Technology, Yokohama, Japan Find articles by Hideki Kobayashi 5 , Dennis Baldocchi Dennis Baldocchi 6 Department of Environmental Science, Policy and Management, University of California, Berkeley, California, USA Find articles by Dennis Baldocchi 6 Author information Article notes Copyright and License information 1 Digital Forest Lab, Department of Geomatics Sciences, Laval University, Quebec City, Quebec, Canada 2 NSF National Center for Atmospheric Research, Boulder, Colorado, USA 3 NOAA/Air Resources Laboratory, Boulder, Colorado, USA 4 School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts, USA 5 Institute of Arctic Climate and Environment Research, Japan Agency for Marine‐Earth Science and Technology, Yokohama, Japan 6 Department of Environmental Science, Policy and Management, University of California, Berkeley, California, USA * Correspondence: Martin Béland ( [email protected] ) ✉ Corresponding author. Revised 2026 Mar 13; Received 2025 Nov 26; Accepted 2026 Mar 18; Issue date 2026 Apr. Published 2026. This article is a U.S. Government work and is in the public domain in the USA. Global Change Biology published by John Wiley & Sons Ltd. This is an open access article under the terms of the http://creativecommons.org/licenses/by-nc-nd/4.0/ License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made. PMC Copyright notice PMCID: PMC13090808  PMID: 41999153 ABSTRACT Multilayer canopy models have been developed for several decades, and interest in this model class remains high despite their complexity because of their multiple advantages over big‐leaf models. Nowadays, a number of scientific and technological advancements favour improving and evaluating these models. First, the advent of lidar technology has enabled detailed mapping of leaves and stems in forest canopies and enhanced the modelling of absorbed solar radiation by both elements within a canopy. Second, long‐time series of eddy covariance measurements are available for model evaluation, capturing years with extreme conditions like droughts. Here we present a new version of the CanVeg model (CanVeg2) with three main modifications: (1) radiative transfer modelling using 3D ray tracing and ground lidar, (2) the addition of a stem energy balance module and (3) the use of a higher order closure model to provide profiles of horizontal wind velocity and variance in vertical wind velocity. We evaluate the CanVeg2 model at five broadleaf forest sites with contrasting canopy structures using long term eddy covariance records, and at two of the sites, soil temperature profiles and scalar vertical profiles. The modelled latent heat and CO 2 flux densities closely matched the eddy covariance measurements. The sensible heat flux density had lower coefficients of determination across sites. The modelled soil heat flux densities were notably higher than the measurements at the three sites where this flux is measured. Divergences between modelled and measured scalars in the lower canopy layers suggest the model would benefit from improved information on soil‐litter moisture content. The favourable results obtained in the model evaluation suggest it may be useful towards addressing new science questions by enabling the estimation of certain canopy states nearly impossible to measure at the canopy level like vertically resolved leaf and stem temperatures and to study the interactions between highly interconnected biophysical and physiological processes. Keywords: ecophysiology, ecosystem fluxes, eddy‐covariance measurements, forest canopy structure, multilayer canopy model The CanVeg2 biophysical model is presented with a focus on three novel components with regards to the CanVeg model: 3D ray tracing radiative transfer modeling, a stem energy budget, and wind and turbulence profiles from higher order closure modeling. The model is evaluated against measurements of turbulent and CO 2 fluxes, net radiation, longwave emissions, as well as vertical profiles of air CO 2 , temperature and humidity, and soil temperature. 1. Introduction Forest canopy functioning models are commonly used to estimate heat, water and CO 2 fluxes between a forest and the above atmosphere. Multilayer models are significantly more complex than “big‐leaf” single‐layer models—used in most large spatial scale land surface models—but are considered more accurate (Bonan et al. 2021 ). Science questions related to interactions between canopy microclimate and leaf physiology require the details on within‐canopy processes provided by multilayer canopy models (Raupach and Finnigan 1988 ). Further, there are instances where big‐leaf and multilayer model estimates of whole canopy level fluxes differ significantly, largely because non‐linear responses to environmental forcings are not considered when treating the canopy as a bulk volume (Beyschlag and Ryel 2007 ; Rastetter et al. 1992 ). Multilayer canopy models have a long development history going back more than half a century, and, besides better correctness, they present advantages over big‐leaf models like permitting the prediction of temperature and humidity profiles within a canopy (Cowan 1965 ; Waggoner and Reifsnyder 1968 ). There are thus situations where one may opt for a multilayer over a big‐leaf model. One such situation is when the effects of non‐linearity in the equations representing the physical and physiological processes on the whole canopy fluxes are important in their own right (Raupach and Finnigan 1988 ). For example, non‐linearity in the light environment can be captured with a sunlit/shaded big‐leaf canopy (dePury and Farquhar 1997 ; Wang and Leuning 1998 ). However, the response of a canopy to extreme events such as heat and drought, which manifests through vertical variation in leaf temperature, leaf water potential and stomatal conductance, cannot be represented in a one‐layer big‐leaf canopy without resorting to empirical prescriptions. In a broad sense, multilayer canopy models hold potential for (1) improving ecosystem response predictions in large scale land surface models, (2) serving as a tool at specific sites to estimate canopy level variables which are impractical to measure and (3) testing our understanding of physical‐physiological processes and their interactions by comparing model estimates against eddy covariance flux tower observations. Multilayer forest canopy models have largely evolved since the 1990s from developments in crop micrometeorological and physiological models in the 60s through the 80s (Goudriaan 1977 ; Jarvis et al. 1985 ). The adaptation to forests required a characterisation of canopy structure and its influence on the radiative transfer, which motivated development of the MAESTRO/MAESTRA model (Medlyn 2004 ; Wang and Jarvis 1990 ). Canopy structure also influences turbulent transport and scalar profiles within the canopy (Baldocchi and Meyers 1988 ; M. Raupach 1988 ), and this body of work led to models like CUPID (Norman 1979 ; Wilson et al. 2003 ) and CANOAK/CanVeg (Baldocchi and Harley 1995 ; Harley and Baldocchi 1995 ). None of these models explicitly considered the role of wood structures in the absorption of solar radiation and heat storage. Recently, the amount of radiation absorbed by stems has been quantified (Béland 2025 ), and the role of this absorption on heat storage and air temperatures diurnal patterns has been suggested to be significant (Swenson et al. 2019 ). Béland ( 2025 ) showed that stems account for about 30%–35% of the near infrared radiation (NIR) absorbed within deciduous forests, indicating that adding stems energy balance to multilayer models may improve canopy microclimate simulations. However, few land surface models explicitly consider the role of wood structures in radiation absorption and heat storage, and in its current version, the Community Land Model (CLM) uses a simple ad‐hoc approach to partitioning absorbed radiation by leaves and wood (Swenson et al. 2019 ). Forests often have very complex structures, with leaves clumped at different scales and in specific areas of the canopy (Béland and Baldocchi 2020 ) and wood structures that widely vary in size in the canopy space. Recent use of ground lidar improved our capability to derive 3D information on canopy structure and differentiate leaves from wood structures therein. Such lidar‐derived forest digital replicas provide inputs for ray tracing radiative transfer models to improve estimates of radiation fluxes absorbed by leaves and wood in 3D, accounting for both the vertical and horizontal heterogeneity of forest canopies using voxel arrays (Béland and Kobayashi 2024 ). Kobayashi et al. ( 2012 ) used the multilayer canopy model CanVeg with radiative transfer computed from geometrical shapes (ellipsoids) derived from airborne lidar in a savanna, but ground lidar and voxel arrays have not yet been used to compute radiative fluxes within a multilayer canopy functioning model. Better data on canopy structure derived from lidar can also improve the computation of vertical variations in wind drag from leaves with higher order closure models to estimate wind velocity and variance in vertical velocity. These models can reveal secondary wind maxima in the stem space of a forest, which can produce more representative dispersion matrices that are used to compute scalar fields and how those scalar fields feed back onto the vertical profiles of scalar and energy fluxes. Improved vertical profiles of vertical wind velocity variance are needed to provide better information for the Lagrangian dispersion matrices that are in turn needed to assess the effects of non‐local turbulent transport in canopies, that is insufficiently represented from K‐theory (Raupach and Finnigan 1988 ). In this paper, we present modifications to the CanVeg multilayer model of Baldocchi and Harley ( 1995 ) to (1) integrate radiation fluxes on leaves and stems computed in 3D using the FLiESvox radiative transfer model (Kobayashi and Iwabuchi 2008 ), (2) compute stem temperatures, heat storage and longwave emissions and (3) compute vertical profiles of wind velocity and vertical turbulent mixing from the Higher Order Closure model of Meyers and Tha Paw U ( 1986 ) and integrate the turbulent term into CanVeg's Lagrangian dispersion matrix. We ran the modified model—hereafter referred to as CanVeg2—using canopy structure derived from ground lidar collected at three broadleaf deciduous forest sites. Each plot surveyed is located within the footprint of a flux tower. We used meteorological records from the towers as forcing to the CanVeg2 model, and the tower measured fluxes over multiple years to evaluate the modelled energy and water fluxes. We also used tower‐based measurements of outgoing longwave radiation flux, air temperature, CO 2 and relative humidity profiles, as well as soil temperature measured at different depths as validation points for the model. 2. Materials and Methods 2.1. The CanVeg Model Description CanVeg is a multilayer canopy functioning model developed to estimate photosynthesis and evaporation rates, as well as sensible heat flux. It has also been used to compute stable carbon isotopes (Baldocchi and Bowling 2003 ) and isoprene emissions (Baldocchi et al. 1999 ), and to study the effect of diffuse radiation on forest canopy fluxes (Knohl and Baldocchi 2008 ), the relation between photosynthesis and radiation reflected by a crop canopy (Baldocchi et al. 2020 ) and the relation between photosynthesis and forest canopy temperature (Helliker et al. 2018 ). The model was initially presented by Baldocchi and Harley ( 1995 ) and Harley and Baldocchi ( 1995 ). The CanVeg model is not coupled to an atmospheric model, it is forced using meteorological measurements usually made at the top of eddy covariance flux towers. Those measurements are recorded every hour or half‐hour, and for forests include incoming shortwave radiation (W m −2 ), air temperature (K), vapor pressure deficit (Pa), wind velocity (m/s), CO 2 concentration (ppm), atmospheric pressure (kPa), friction velocity (m/s), soil temperature (K) and soil moisture (fraction). The model uses vertical profiles of leaf area index (LAI), foliage clumping and leaf angle distribution function to describe the canopy structure. It also uses numerous leaf level physiological parameters as well as optical properties for leaves and soil, these parameters and values used will be presented in Section 2.2 . The CanVeg model has three main components: biophysical, physiological and turbulent transport. The biophysical component computes the radiative transfer of solar radiation through the canopy layers to estimate the amount of photosynthetically active radiation (PAR) and near infrared radiation (NIR) absorbed by sunlit and shaded leaves in each layer and the leaf and soil energy balance which leads to estimates of sunlit and shaded leaf and soil temperatures. A radiative transfer process is also used to model the emission and absorption of longwave radiation based on the leaf and soil temperatures. The physiological component computes the leaf stomatal conductance, the photosynthesis and transpiration rates and sensible heat flux, which are non‐linear functions of the biophysical variables. The turbulent transport component computes the vertical profiles of CO 2 concentration, air temperature and vapor pressure based on the sources and sinks of CO 2 , sensible heat and latent heat in each canopy layer using a dispersion matrix developed from Lagrangian theory that releases an ensemble of particles into the air space (M. Raupach 1988 ). Several of the processes involved in these three above components are highly coupled. For example, leaf stomatal conductance, photosynthesis rate, leaf temperature and leaf transpiration are all interdependent. The initial CanVeg model was revolutionary at its creation in resolving these interdependencies by iterating through part of the calculations until leaf temperatures stabilise between successive iterations. In the first iteration, the longwave radiative transfer is first performed assuming leaf temperatures equal air temperature. The leaf and soil energy balance then adjusts leaf and soil temperatures based on the radiative forcings and the turbulent diffusion of heat. First estimates of stomatal conductance, photosynthesis and transpiration rates and heat exchange from leaves are made and the turbulent transport computes vertical profiles for the CO 2 , air temperature and vapor pressure scalars. The subsequent iterations start with longwave radiative transfer using the updated leaf and soil temperatures and follow through to the turbulent transfer. Five to ten iterations are usually sufficient for the average of all leaf temperatures to stabilise. In the following subsections we will describe the mathematical treatments used to express each process in the CanVeg2 model. The CanVeg2 model code retains the same model structure as the CanVeg model presented in Baldocchi and Harley ( 1995 ), while several of the code functions were replaced with adapted functions developed by Bonan ( 2019 ) and used in the CLM‐ml multilayer model (Bonan et al. 2021 , 2018 ). These functions include the 1D radiative transfer modelling for longwave radiation, the leaf energy balance, plant hydraulics and stomatal conductance, photosynthesis and the soil energy balance functions. The turbulent transport function using a Lagrangian dispersion matrix presented in Baldocchi and Harley ( 1995 ) is maintained. The fundamental advancements between CanVeg and CanVeg2 relate to (1) the radiative transfer process being done in 3D from ground lidar data and including stem radiation absorption, (2) the stems being considered in the canopy energy balance and emitting longwave radiation and (3) the wind velocity and turbulence statistic profiles are derived using the Higher Order Closure model of Meyers and Tha Paw U ( 1986 ). The higher order closure principles are used here only to compute the mean horizontal wind profiles and variance in vertical wind speed profiles as a function of LAI profiles that drive the Lagrangian dispersion matrix. We chose not to use a complete higher order closure model to simulate turbulent diffusion, which would require modelling conservation budgets and resolving second and third order terms as done by Meyers and Tha Paw U ( 1987 ) and Pyles et al. ( 2003 ). The use of the higher order closure model is further described in Sections 2.1.1 and 2.1.7 below. A multilayer canopy computer model involves representing processes with equations, correctly representing the dependencies between these processes and developing numerical methods to solve the equations in a way that respects the interactions and feedbacks between processes. In the next sections we attempt to present the main equations involved and the numerical methods used to solve them. The order in which we present the processes being modelled follows, for the most part, the order in which they are treated in the model. 2.1.1. Wind Velocity and Turbulence Vertical Profiles Vertical profiles of wind velocity are a critical component of the model, since they influence the turbulent diffusion of heat and thus significantly influence leaf temperatures, a central variable to the model considering its relation to multiple processes. The wind velocity and the variance in vertical wind velocity profiles were computed for each site using the Higher Order Closure model of Meyers and Tha Paw U ( 1986 ). The profiles are computed for neutral conditions on the basis of canopy height, LAI vertical profile and canopy drag coefficient, which encompasses spatial scales from leaf to canopy. Using this new information on canopy vertical profiles of LAI from lidar measurements in higher order closure models gives us better information on the turbulence fields of wind velocity and variance in vertical velocity. This information, in turn, is needed to refine the computations of dispersion matrices that are used to compute profiles of scalars that both drive local fluxes and are the result of those fluxes. The two vertical profiles are produced once for each site, read by the model which then scales the profiles for different thermal stability conditions (further details are provided in Section 2.1.7 below). Baldocchi and Meyers ( 1988 ) showed a very close correspondence between the measured mean wind velocity profile in a deciduous forest near Oak Ridge, TN, and the Meyers and Tha Paw U ( 1986 ) model. 2.1.2. Radiative Transfer Modelling The partitioning of incoming shortwave radiation between direct and diffuse light is done following Weiss and Norman ( 1985 ), which offers reasonable performance even though better more recent models exist (Oliphant and Stoy 2018 ). The incoming longwave radiation is calculated following Choi et al. ( 2008 ). The shortwave radiative transfer within the canopy can be simulated using either the 1D model from Norman ( 1979 ) (modified by Béland et al. in review to include stems' radiation absorption), or the 3D FLiESvox ray tracing model presented in Kobayashi and Iwabuchi ( 2008 ) and Béland and Kobayashi ( 2024 ). The latter relies on voxel‐based 3D canopy structure information derived from ground lidar (methods described in Béland et al. ( 2011 ) and Béland et al. ( 2014 )), where a given voxel within the array may contain different densities of leaves, wood or both elements. The interception of radiation by stems using this model has been studied in Béland ( 2025 ). Both 1D and 3D models consider vertical profiles of foliage clumping and leaf angle distributions. The longwave radiative transfer is done in 1D using the Norman ( 1979 ) model modified by Béland et al. ( in review ) to include radiation interception and emission by stems. The 3D radiative transfer model simulates the interaction of photons with leaves, stems and the soil by tracing the path of a large number of parcels emitted from a direct (sun position) or diffuse source. The fate of a given parcel upon interaction depends on the optical properties of the intercepting element; it can be reflected, transmitted or absorbed. The FLiESvox model is too computationally demanding to be run at each time step of the CanVeg2 model; hence, we used a look‐up table approach where FLiESvox is run at each site using nine different sun zenith angles (5° to 95° with 10° intervals) and nine diffuse light fractions for a total of 81 model runs for PAR and 81 runs for NIR. For each run, the absorbed radiation fluxes calculated at the voxel level for leaves and wood are summed horizontally to yield 1D vertical profiles of absorbed radiation. The absorbed radiation profiles corresponding to the closest sun zenith angle and diffuse light fraction are then read in as input into CanVeg2's leaf and stem energy balance functions described in the next sections. 2.1.3. Leaf Energy Balance, Stomatal Conductance, Photosynthesis, Transpiration and Respiration The water‐use efficiency optimisation theory formulated by Cowan ( 1978 ) is used to calculate the stomatal conductance which maximises the marginal carbon gain (through photosynthesis) of water loss (through transpiration) over the model time step, which corresponds to the interval over which the eddy covariance measurements are averaged and recorded (typically 30 min or 1 h). A site‐specific parameter called the marginal water‐use efficiency (𝜄) is used to relate transpiration cost to a carbon gain, effectively setting the relative cost of transpiration: ∂ A n ∂ g sw ∂ E ∂ g sw = ι (1) where A n is the CO 2 assimilation rate, E is the transpiration rate and g sw is the stomatal conductance. Here we assume that 𝜄 is constant through time, meaning this cost–benefit analysis is unaffected by changing conditions like air CO 2 concentration and/or VPD throughout the day. The iterative numerical methods presented in Bonan et al. ( 2014 ) for CLM‐ml (Bonan et al. 2021 , 2018 ) are used to solve the system of equations related to leaf temperature, transpiration and photosynthesis so that further stomatal opening yields insufficient carbon gain per water loss. The leaf photosynthesis model is that of Farquhar et al. ( 1980 ), with the co‐limitation between Rubisco limited and RuBP regeneration of Collatz et al. ( 1990 ) to smooth the transition. The leaf physiology parameters used are listed in Supporting Information , Table S1 . Rates of maximum carboxylation, electron transport and respiration were corrected for their dependence on leaf temperature. Sensible and latent heat fluxes require calculation of the leaf boundary layer conductance to heat and water respectively. The boundary layer conductances are calculated by summing the forced and free convection regimes, thus assuming they occur together, and follow the description provided in Bonan ( 2019 ), with the added effect of leaf clumping on the laminar and turbulent air flows. Landsberg and Powell ( 1973 ) showed that when tree leaves are clumped, the mutual interference significantly increases the boundary layer resistance. Baldocchi ( 1991 ) pointed to the work of (Grace and Wilson 1976 ) to suggest that the Sherwood and Nusselt numbers for heat and water vapor may be twice the number derived from flat plate theory, thus decreasing boundary layer resistance, and that the combined effect with that of mutual interference may cancel out in forest canopies. Since we are considering here vertically heterogeneous profiles for leaf clumping, we are adjusting the Sherwood and Nusselt numbers so that the boundary layer resistance is higher at the canopy tops where leaves are very clumped. This is done by multiplying both numbers by the leaf clumping factor + 0.35 for a given canopy layer. The clumping factor in the lower canopy layers is about 0.65, and the clumping of leaves at the branch scale was shown to be absent in those layers by Béland and Baldocchi ( 2021 ) (clumping occurring at larger scales is assumed not to influence leaf boundary layer resistance), hence the multiplicative factor for the effect of mutual interference is one. At canopy tops where leaf clumping factors are about 0.35 in temperate broadleaf forests, the multiplicative factor applied to the Sherwood and Nusselt numbers is 0.7 which increases boundary layer resistance. For all canopy layers, the Sherwood and Nusselt numbers are multiplied by an empirical correction factor of 1.5 suggested by Schuepp ( 1993 ) to account for leaves having lower resistances than a flat rectangular plate. We imposed a minimum value of 0.2 mol m −2 s −1 for leaf boundary conductance to heat to avoid unrealistic conditions in very low wind speeds (below 0.1 m s −1 ), similarly to Bonan et al. ( 2026 ). The sensible and latent heat flux calculations follow the mathematical development of chapter 10 in Bonan ( 2019 ). The leaf temperature is the value at which the leaf energy budget is balanced. Terms considered in the budget are radiative forcing (shortwave and longwave radiation absorbed less emitted longwave radiation), sensible and latent heat, and the leaf heat storage (which is very small relative to other terms). Sunlit and shaded leaf temperatures are calculated by solving the energy budget using the Newton–Raphson iteration method presented in Bonan ( 2019 ) chapter 10. Leaf respiration is calculated as a function of leaf temperature using a peaked Arrhenius function as described in Bonan ( 2019 ) chapter 11. Leaf respiration at 25°C was set at 1.5% of Vcmax at 25°C, which value is site‐specific (see Table S2 in Supporting Information ). 2.1.4. Stem Energy Balance The stem surface temperature and the temperatures of multiple layers within the stems are calculated using an energy balance approach, where the net radiation is balanced by sensible heat, heat conduction and heat storage. The heat flux into the stem is: F 0 T = Q a − εσ T 4 − c p T − T air g ac (2) where T is the stem surface temperature, Q a is the radiative forcing on stems (the sum of absorbed PAR, NIR and longwave radiation), the second term on the right is the emitted longwave radiation, and the third term is sensible heat, with g ac the resistance to convective heat transfer. In the 3D radiative transfer model, wood area is treated as a turbid medium within the voxel cubes. Hence, the wood structures do not have an azimuthal orientation and no sunlit or shaded side. The wood structure is treated as a cylinder only in the simulation of the convective and conductive heat exchanges. The calculation of heat loss from stems by moving air uses a resistance to convective heat transfer over a cylindrical object exposed to cross‐diameter air flow following Monteith and Unsworth ( 2013 ). For forced and free convection, the Reynolds and Grashof numbers used to calculate the Nusselt number are influenced by the stem diameter, and their calculation follows table A5 in Monteith and Unsworth ( 2013 ). Nusselt numbers for forced and free convection are summed, assuming they occur in parallel. The stem boundary layer conductance for heat (mol m −2 s −1 ) is obtained by multiplying the Nusselt number by the molecular diffusivity of heat divided by the stem diameter at the given canopy layer. The Nusselt number is multiplied by an adjustment factor of 0.5 to account for the fact that stems have a non‐flat surface, with bark roughness influencing air flow. A thin first layer 0.0001 m in thickness is used to represent the stem surface conditions. The stem cylinder volume is further divided in 10 layers, each having equal volume (pi*h*R 2 /10, where h is the cylinder length corresponding to the canopy layers thickness (30 cm) and R is the cylinder radius, see Figure 1 ). The heat storage in the stem layers is estimated by modelling the heat transfer within the boles following radiative forcing and surface heat exchange. The stem surface heat exchange is solved simultaneously with within stem temperatures using the implicit method presented in Bonan ( 2019 ) section 7.3. This method uses a Taylor series approximation to solve the stem surface temperature at time t + 1 using the fluxes and surface temperature at time t following: F 0 T n + 1 = F 0 T n + ∂ F 0 T n ∂ T T n + 1 − T n (3) where F 0 is the net energy flux into the stem (W/m 2 ), T is the stem surface temperature and n is time. The difference between time n and n + 1 corresponds to the time elapsed between measurement records at the tower, which is also the model time step, either 30 min or 1 h. The longwave and sensible heat fluxes and the heat storage are thus first calculated using stem surface temperature at time t ; the within stem layers temperatures are updated based on the energy flux into the stem and the conduction of heat between layers, which leads to an updated stem surface temperature. The longwave and sensible heat fluxes and the heat storage are then calculated using the stem surface temperature at time t + 1. FIGURE 1. Open in a new tab Illustrations of a stem cross section with 10 layers used to calculate heat conductance within the stems, where each layer has the same area (or volume if the full cylinder is considered). In the first moment of the first model iteration, the temperatures of all stem layers are initialised at the mean daily air temperature of the month being processed and thereafter are calculated from the stem energy budget. The centre of the stems is a boundary condition with zero heat flux. The first thin stem layer is the other boundary condition to calculate the heat transfer using Fourier's law: c v ∂ T ∂ t = ∂ ∂ r κ ∂ T ∂ r (4) where c v is the heat capacity, 𝜅 is the thermal conductivity, ∂ T ∂ r is the radial temperature gradient and ∂ T ∂ t is the change in temperature with time. The numerical solution to solve the temperature of each stem layer uses a tridiagonal system of equations presented in Bonan ( 2019 ), chapter 5. The stem emissivity and thermal properties used are provided in Supporting Information , Table S1 . Since the amount of heat a stem can store is a function of stem diameter, we use a mean stem diameter vertical profile defined from ground lidar data. The stem sensible heat flux and heat storage flux are thus calculated using the average stem diameter for a given canopy vertical layer with units W m −2 of wood area, which is then converted to W m −2 of ground area by multiplying by the layer wood area index. Because larger trees represent a disproportionally large fraction of the total wood biomass in a forest, only the trees reaching a height above about two thirds of the maximum canopy height are considered in the calculation of the average stem diameter profile. To determine the effect of using an average stem diameter on the stem heat storage calculations, we ran a model simulation at the EMS site calculating the heat storage from three stem diameter profiles corresponding to three classes of tree heights: (1) trees having between 5 and 15 m height, (2) 15–21 m in height and (3) 21–25 m in height. We then calculated the whole canopy stem heat storage in W m −2 of ground area using the wood area index of each class and compared it to the storage flux calculated using the average stem diameter profile. The differences between the methods were minimal which suggested it is suitable to use a single stem diameter profile representing the mean stem diameter per vertical layer of the larger trees in the canopy. Further details on the diameter profiles used are provided in Section 2.2.1 below. 2.1.5. Plant Hydraulics Plants may adjust stomatal conductance to reduce water loss, the risk of hydraulic failure and loss of conductance following a decrease in xylem water potential due to low soil matric potential and plant hydraulic architecture. In CanVeg2, the leaf water potential is calculated for each vertical layer and each model time step following: ψ = ψ s − ρ wat ⋅ g ⋅ h − E K L (5) where ψ s is the soil water potential (Pa), ρ wat ·g·h is the gravitational potential at height h above ground, E is the transpiration flux and K L is the whole‐plant hydraulic conductance (mol H 2 O m −2 s −1 Pa −1 ) which integrates the hydraulic conductance of roots, stems and branches. K L was calculated following the same approach and parameter values as in Bonan et al. ( 2018 ). The stomatal conductance may be reduced if the leaf water potential falls below a hydraulic safety threshold. This threshold, as well as the amount by which stomatal conductance is reduced for a given value of leaf water potential depends on the mix of hydraulic strategies of tree species within a given plot being modelled, that is, either isohydric (more conservative strategy where stomatal conductance is reduced more rapidly upon water stress) or anisohydric (riskier strategy where high stomatal conductance is maintained as soils dry and/or vapor pressure deficit increases). Following the calculation of the optimal stomatal conductance described in Section 2.1.3 , we use an approach which allows a smooth and gradual reduction of stomatal conductance caused by the hydraulic failure risk constraint if leaf water potential falls below the hydraulic safety threshold. A single two‐parameter function presented in Christoffersen et al. ( 2016 ) is used to represent the decline in stomatal conductance in relation to leaf water potential for an average strategy representative of the ensemble of species within the flux tower footprint. The equation is. g sw = g sw 1 + ψ ψ 50 a − 1 (6) where ψ is the leaf water potential, ψ 50 is the leaf water potential at 50% loss of conductivity, a is a shape parameter. The optimality model for stomatal conductance presented in Section 2.1.3 therefore has two criteria determining the conductance value, the water‐use efficiency and hydraulic safety. This approach is inspired by the Soil–Plant‐Atmosphere model of Williams et al. ( 1996 ) presented in Bonan ( 2019 ), chapters 12 and 13 and used in CLM‐ml (Bonan et al. 2021 , 2018 ). Leaf water potential is largely influenced by the soil matric potential, and any attempt to simulate the response of stomatal control to low soil water content must rely on an accurate description of the relation between soil water content (which is measured at flux tower sites) and soil matric potential. This relationship is highly influenced by soil type and can vary drastically between sites. This will be further discussed in the next section and Section 2.2.2 below. 2.1.6. Soil Energy Balance The soil is represented using ten vertical layers of different thicknesses. All layers are assumed to have the same water content, corresponding to the field measured value generally made at a depth of about 15 cm. Even though the amount of precipitation is generally recorded at flux tower sites, we did not include modelling of the infiltration of rainwater into the soil. The soil matric potential is calculated from the measured soil moisture content using an equation derived from pre‐dawn leaf water potential measurements at sites where they are available. The temperatures of soil layers are initialised using the tower recorded soil temperature (measured at about 15 cm depth) for layers down to 40 cm depth, and the mean yearly above canopy air temperature for deeper layers between 40 cm and 3 m. In subsequent iterations the measured and modelled soil temperature are allowed to diverge, enabling the use of the measured soil temperature as a validation point to assess the modelled values. An energy balance approach with multilayer heat transfer by conduction, similar to the one used for stem energy balance, is used to calculate soil temperatures. The radiative forcing is balanced by turbulent heat exchanges between the soil and the air above, and by the conduction of heat within the different soil layers. The Taylor series approximation presented in Bonan ( 2019 ), chapter 7, is used to solve for soil surface temperature. The soil evaporation is modelled using the approach presented in Baldocchi et al. ( 2000 ), where the soil surface aerodynamic resistance accounts for thermal stratification and the calculation follows the method presented by Daamen and Simmonds ( 1996 ), and the soil resistance to evaporation used the empirical relationship from Sellers et al. ( 1996 ). In defining the soil thermal properties, the ten vertical layers are separated in three strata: litter (between 0 and 2 cm), organic (between 2 and 18 cm) and mineral (deeper than 18 cm). All soil layers' thermal conductivity and heat capacity were calculated as an average of solid, air and water fractions for each stratum following Bonan ( 2019 ) chapter 5. 2.1.7. Turbulent Transport The turbulent dispersion of the calculated scalar atmospheric constituents heat, water vapour and CO 2 is modelled using a Lagrangian dispersion matrix based on M. R. Raupach ( 1989 ) and as implemented by Baldocchi and Harley ( 1995 ). Lagrangian models simulate the transport process by tracking the trajectories of an ensemble of fluid parcels as they are advected and diffused by the mean wind and turbulence (Baldocchi and Meyers 1989 ). The dispersion matrix uses a vertical profile of the standard deviation of the vertical wind velocity (𝜎 w ), which is calculated using the Higher Order Closure model of Meyers and Tha Paw U ( 1986 ). The scalar source strengths and concentration profiles are related by: c i − c ref = ∑ j = 1 N S c , j Δ z j D ij (7) where c i is the mean scalar concentration, c ref the scalar concentration at a reference level (i.e., the tower measurement height), N is the canopy layer, S c,j 𝛥 z j is the source flux at level z and D ij is the dispersion matrix. A look up table approach is used to account for different atmospheric stability conditions, where seven different dispersion matrices ( D ij ) are produced for different Obukhov length scale values corresponding to stable, unstable and neutral conditions. At each model time step, the Obukhov length scale is calculated from the vertically integrated sensible heat fluxes and friction velocity. Each stability condition is associated with a 𝜎 w value above the canopy, and the entire 𝜎 w profile modelled for neutral conditions is scaled to fit the above canopy value for the various stable and unstable conditions. 2.2. Model Evaluation As stated by Medlyn et al. ( 2005 ), it is useful when evaluating a model to define the purpose or purposes the model is intended for. Here we consider that getting the correct fluxes for the correct reasons may be at the forefront of concerns when evaluating the capacity of the model to match observed fluxes. We opted to keep to a minimum the number of parameters adjusted at each site (four physiological parameters related to photosynthesis and stomatal conductance), and those parameters not to vary with time. Using the same model parameters across different years and environmental conditions (including extremes) may indicate that the model is getting the right answer for the right reasons when validating against eddy covariance data. Another intended purpose of the model is the prediction of canopy response to extreme conditions, like high heat combined with drought. This requires, among other things, reliable simulation of difficult‐to‐measure quantities like leaf temperatures and their spatial variability. We thus included additional validation checks using outgoing longwave radiation flux measurements made from the flux towers and soil temperature profiles available at one of the sites used. It also requires a correct representation of stomatal response to low soil water content and high vapor pressure deficit, as well as reasonable wind velocity profiles that drive the turbulent heat exchanges. The CanVeg2 model is forced using meteorological data recorded at the flux towers every half‐hour or 1 h. Gap‐filled data records were used at all sites. Those measurements made above canopy are: air temperature, incoming shortwave radiation, vapor pressure deficit, wind speed, air CO 2 , atmospheric pressure, friction velocity and below canopy at a depth of about 15 cm: soil temperature and soil water content. The model outputs CO 2 fluxes, and sensible and latent heat fluxes. At one site (the Harvard Forest Environmental Monitoring Site, described below), vertical profile measurements of air temperature and relative humidity were compared against the model estimates to evaluate the accuracy of the sources/sinks of heat and water vapor and their vertical mixing through the dispersion matrix. 2.2.1. Sites Three sites are used in this study, representing a gradient in canopy structure from a smooth canopy top in a temperate forest to a very complex old growth tropical forest. We also included one site at which a drought/heat event occurred during the period covered by flux measurements. Each site was surveyed with ground lidar over a plot located within the flux tower footprint to derive accurate descriptions of canopy structure in 3D space, including leaf and wood area. Methods used to acquire and process the lidar data are described in detail in Béland et al. ( 2014 ), Béland and Kobayashi ( 2021 ) and Béland ( 2025 ). The leaf area density is estimated at the voxel level (30 cm in size length) based on the interception and transmission of laser pulses through each voxel volume. Foliage clumping at the within voxel scale was estimated using vertical profiles following the relation with cumulative leaf area index presented in Béland and Baldocchi ( 2021 ). The leaf area index (LAI) and foliage clumping factor profiles for each site are shown in Figure 2 . The stem mean radius profiles used in the calculation of stem heat storage mentioned in Section 2.1.4 above are shown in Figure 3 . FIGURE 2. Open in a new tab Canopy structure represented using voxel cubes (left), with colors representing height above ground and leaf area index (middle) and within voxel foliage clumping profiles (right). This figure translates that foliage clumping is considered at two spatial scales in the radiative transfer function: first at the branch scale (right) and at large spatial scales associated with the relative positions of tree branches and crowns with regards to a given sun illumination direction (left). However, when considering the effect of foliage clumping on leaf boundary layer conductance, only the branch level clumping (right) is relevant, since it is the clumping of leaves in close proximity that influences the air flow on the leaf surfaces. FIGURE 3. Open in a new tab Mean stem radius profiles for all sites used in the calculation of stem heat storage. The first site is at the Harvard Forest Environmental Monitoring Site (EMS), Petersham, Massachusetts, USA (42° 32′ N, 72° 10′ W, elevation 340 m). A ground lidar survey was performed on a 60 × 60 m 2 plot within the flux tower footprint in August 2017. The average canopy height is about 25 m with a LAI estimated at 4.8. The tower measurement height is 30 m. The selected plot is dominated by red oak ( Quercus rubra ) and red maple ( Acer rubrum ). The leaf angle distribution is planophile up to the canopy top (Béland and Kobayashi 2021 ). The forest stand is about 110 years old, and most of the foliage is concentrated towards the top (top‐heavy canopy), with a very flat canopy surface. The mean yearly air temperature is about 9°C and the mean annual precipitation is 1071 mm. The soil texture used in CanVeg2 for this site is loam. The eddy covariance flux tower at this site records meteorological and flux data hourly since 1992 (Munger 2022 ). The CanVeg2 model was run for the month of July at this site between 1993 and 2020. Several additional available measurements at this site were used to evaluate the scalar vertical profiles produced by the model (Munger and Wofsy 2024 ). First, measurements of CO 2 concentration, air temperature and relative humidity were made from different heights above ground. The CO 2 concentrations were measured at heights 0.3, 0.8, 4.5, 7.5, 12.7, 18.3, 24.1 and 29 m between 1993 and 2004. Relative humidity and air temperature were measured at heights 2.5, 7.6, 15.4, 22.6 and 27.9 m between 2012 and 2020. Second, soil temperature measurements were made between 2015 and 2020 from three locations near the plot used. At each location soil temperature was recorded at three depths: 5, 20 and 50 cm. Recordings from the three locations were averaged for each depth to provide spatially representative soil temperature profiles. The second site is in the Morgan‐Monroe State Forest, near Bloomington, Indiana, USA (39°19′ N, 86°25′ W, elevation 250 m). A ground lidar survey was performed on a 60 × 60 m 2 plot within the flux tower footprint in July 2018. The average canopy height is about 30 m and the LAI is about 5.1. The tower measurement height is 46.2 m. The plot is 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 upper canopy layers (Béland and Kobayashi 2021 ). The soil texture used in CanVeg2 for this site is loam. The mean yearly air temperature is about 12.5°C and the mean annual precipitation is 1032 mm. The site has an eddy covariance flux tower that is part of the Ameriflux network and records hourly data since 1999 (Novick and Phillips 2022 ). This site had a significant heat and drought event in 2012, with the peak of temperatures occurring in July. The CanVeg2 model was run for the month of July between 2000 and 2020. The third site is in the Pasoh Forest Reserve of the Forest Research Institute Malaysia, near Simpang Pertang, Malaysia (2°58′ N, 102°18′ E, elevation 140 m). A ground lidar survey was performed on a 60 × 60 m 2 plot within the flux tower footprint in April 2018. The average canopy height is 35 m, and the tower measurement height is 54 m. The LAI is about 7.3 (Kira et al. 2013 ), and the canopy is 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 (Kosugi and Takanashi 2019 ). The mean yearly air temperature is about 25°C, and the mean annual precipitation is about 2000 mm. The soil texture used in CanVeg2 for this site is silty clay. This site has an eddy covariance flux tower recording half‐hour data, with data submitted to FLUXNET (Pastorello et al. 2020 ) between 2003 and 2009. It is important to note that the turbulent heat fluxes available for this site in the FLUXNET database were adjusted to close the energy budget using the Bowen ratio method. Hence, the fluxes used in this study for this site are the Bowen ratio corrected values and not the raw measurements. The CanVeg2 model was run at this site for the month of February between 2003 and 2009. The fourth site is the Bartlett Experimental Forest, within the White Mountains National Forest in north‐central New Hampshire, USA (44.0646′ N, 71.2881′ W). A ground lidar survey was performed on a 48 × 48 m 2 plot within the flux tower footprint in August 2023. The maximum canopy height is about 25 m, and the tower measurement height is 26.5 m. The LAI is about 5.3 m 2 /m 2 , and the dominant species in the plot are Red Maple, Yellow Birch, White Ash and American Beech. The leaf angle distributions for this site were estimated using digital photography following the method described in Ryu et al. ( 2010 ), and fits the uniform distribution above 13 m and the planophile distribution below 13 m (Hastings 2025 ). The soil texture used in CanVeg2 for this site is loam. The mean yearly air temperature is about 5.6°C and the mean annual precipitation is 1245 mm. The site has an eddy covariance flux tower part of the Ameriflux network and records half‐hourly data since 2004 (Ouimette et al. 2018 ). The fifth site is the Hubbard Brook Experimental Forest, located within the White Mountain National Forest in central New Hampshire (43.9397′ N, 71.7181′ W). A ground lidar survey was performed on a 48 × 48 m 2 plot within the flux tower footprint in August 2023. The maximum canopy height is about 27 m, and the tower measurement height is 30 m. The LAI is about 5.2 m 2 /m 2 , and the dominant species in the plot are Sugar Maple, Yellow Birch, White Ash and American Beech. The leaf angle distributions for this site were estimated using digital photography following the method described in Ryu et al. ( 2010 ), and fits the uniform distribution above 23 m and the planophile distribution below 23 m (Hastings 2025 ). The soil texture used in CanVeg2 for this site is sandy loam. The mean yearly air temperature is about 6°C and the mean annual precipitation is 1400 mm. The site has an eddy covariance flux tower part of the Ameriflux network and records half‐hourly data since 2016 (Kelsey and Green 2025 ). Data at this site is provided as an AmeriFlux BASE data product, and is not generated using the ONEFlux processing pipeline which filters fluxes and u*. At this site, similarly to the EMS site, vertical profiles of measured scalars are available to evaluate the vertical profiles produced by the model. Measurements of air temperature were made from 1.5, 3, 6, 9, 12, 15, 18, 21, 24, 27 and 30 m above ground. Relative humidity was measured at heights 1.5 and 33 m. Soil temperature was recorded at four depths: 5, 20, 30 and 50 cm. The mean horizontal wind velocity and variance in vertical wind velocity profiles derived from the Higher Order Closure model from Meyers and Tha Paw U ( 1986 ) at each site is presented in Figure 4 . A value of 0.25 was used at all sites for aerodynamic drag coefficient. The same leaf and soil optical properties were used at all sites, based on field measurements made at the Harvard Forest site and values published in Majasalmi and Bright ( 2019 ). These values are listed in Supporting Information , Table S1 . FIGURE 4. Open in a new tab Mean horizontal wind velocity and variance in vertical wind velocity profiles in neutral condition from the Higher Order Closure model for three of the sites used (only three shown for clarity). 2.2.2. Model Parameters The CanVeg2 model has numerous parameters describing aspects of leaves, stems and soil, many of which have very little influence on whole canopy functions, while some have very large influence. To test the generality of the model (in contrast with site‐specific calibration), we allowed only four plant physiology parameters to vary between sites. The first parameter is the maximum leaf carboxylation rate ( V cmax ), which was set to 50 μmol m −2 s −1 at the EMS site, 37 μmol m −2 s −1 at Morgan Monroe, 35 μmol m −2 s −1 at Pasoh, 35 μmol m −2 s −1 at Bartlett and 55 μmol m −2 s −1 at Hubbard Brook. These values were obtained by approximately fitting the diurnal pattern amplitude in tower measured gross primary productivity to the model pattern (see Figure S1 in Supporting Information for an example at the Morgan Monroe site). The second parameter is the marginal water‐use efficiency (𝜄), which had values of 700 at the EMS site, 500 at Morgan Monroe, 800 at Pasoh, 1000 at Bartlett and 700 at Hubbard Brook. These values were roughly adjusted to obtain similar Bowen ratios between flux tower measurements and model estimates (see Figure S2 in Supporting Information for an example at the Morgan Monroe site). The other two parameters relate to the function used to reduce stomatal conductance when leaf water potential triggers a hydraulic safety response (described in Equation 6 ). The EMS and Pasoh sites rarely reach leaf water potentials triggering this function, It is more relevant to dry years at the Morgan Monroe site, which had a value of 8 for parameter a, and value of −2.35 for parameter ψ 50 . These values were determined by adjusting the parameters until a reasonable fit was obtained between measured and modeled latent heat during the 2012 drought/heat year. The values for the other sites are provided in Supporting Information , Table S2 . The values for V cmax , 𝜄, a and ψ 50 are similar to values used in CLM‐ml (Bonan et al. 2021 , 2018 ). The relationships between soil water content and soil matric potential (the soil moisture retention curve) used at the Morgan Monroe site were from the relationship developed specifically at this site by Wayson et al. ( 2006 ). The relations at other sites used generic relations used in the Community Land Model (CLM) (Lawrence et al. 2019 ), which is a function of soil texture. We found it was essential to use the retention curve from Wayson et al. ( 2006 ) at Morgan Monroe since none of the generic CLM relations were enabling a correct response to lowering soil water content. Figure 5 shows that for the relation derived at Morgan Monroe, the matric potential decreases slowly with decreasing soil water content, suggesting that water is loosely held in these soils. The shape of the Wayson et al. ( 2006 ) curve is similar to a sand or loamy sand soil in CLM, but the matric potential is much lower. FIGURE 5. Open in a new tab Relations between soil water content and matric potential for different soil textures used in the Community Land Model (CLM) (calculated from the Campbell ( 1974 ) function) shown with thin lines and with the relation published in Wayson et al. ( 2006 ) for the Morgan Monroe site (MMS). None of the CLM curves represent the relation well enough to be used effectively at Morgan Monroe. Soil respiration was set at a fixed value of 6 μmol CO 2 m −2 s −1 , which approximates the average soil respiration measured at the EMS site at different times of day for 13 days during the month of July in 1998, 1999, 2002 and 2003 (average value was 5.82 μmol CO 2 m −2 s −1 , standard deviation of 1.20 μmol CO 2 m −2 s −1 ; Matthes et al. 2025 ). 2.3. Sensitivity Analysis Since estimating stem heat storage using a multilayer canopy model is a novel approach, we performed a sensitivity analysis to evaluate the relative importance of four model parameters: stem diameter profile, stem boundary layer conductance, stem heat capacity and stem thermal conductance. The first parameter determines the volume of wood within which heat can be conducted. The second drives the convective heat exchanges with the surrounding air at the stem surface. The third and fourth parameters determine the heat flow within the stems and the amount of heat required to change the stem temperature. We performed this analysis at the Morgan Monroe site between 2016 and 2020, using for each parameter half and double the default value. 3. Results and Discussion 3.1. Turbulent Heat Fluxes and Net Radiation Figure 6 shows the scatterplot of measured vs. modelled turbulent heat fluxes at the five sites. The latent heat fluxes have higher R 2 values than for sensible heat fluxes. The modelled latent heat fluxes appear to saturate at high values, particularly at the tropical site, this is likely because evaporation from surfaces (soil‐litter layer, leaves, branches and trunks) following rain events are not currently included in the model. Rainfall amounts are recorded at the sites used, and future versions of CanVeg2 may include a scheme to simulate water interception and surface evaporation. This may also improve modelled soil temperature which is presented in the next section. Heat fluxes are noisier at the Hubbard Brook site, particularly at night (shown in blue), likely the effect of data processing not performed with the ONEFlux pipeline which includes flux and u* filtering. Both turbulent heat fluxes are underestimated by the model for the Pasoh site (Figure 6 bottom). This is explained by the measured fluxes being corrected for energy closure using the Bowen ratio method and are not the raw measured fluxes at this site, as mentioned in Section 2.2.1 . Since the model includes the stem heat storage term, which is not considered when fluxes are corrected using the Bowen ratio, it is expected that the published corrected measurements of turbulent heat fluxes are overestimated at the Pasoh site. FIGURE 6. Open in a new tab Comparison of measured and modeled turbulent heat fluxes, latent heat (LE, left) and sensible heat (H, right) for the three sites used and all years available at each site. Dashed lines are the 1:1 relation, the solid lines are the linear fit. For the Hubbard Brook site data points for which net radiation is less than 50 W/m 2 are shown in blue. The fluxes measured over the course of 1 month for multiple years were aggregated by calculating the average values for each half‐hour or 1 h period over the given month. This averaging is suggested by Baldocchi and Wilson ( 2001 ) to reduce variability resulting from the stochastic nature of processes involved. The diurnal patterns are shown in Figure 7 for the five sites. At the Morgan Monroe site, modelled sensible heat appears overestimated in the morning and underestimated in late afternoon. It is not clear what generated the high levels of sensible heat in the morning. The underestimation of sensible heat in the afternoon could be due to overestimation of latent heat, and since the model is forcing energy closure, sensible heat fluxes are reduced by the high latent heat fluxes. Interestingly, it appears the modelled and measured sensible heat fluxes were in much better agreement at this site between 2011 and 2020 than they were between 2000 and 2010 (see Figure 8 ). The measured sensible heat fluxes are notably lower for the first decade of measurements than for the second. This may suggest that either something changed in the environmental conditions which is not fully considered by the model, or in the instrumentation used, or the measurement and data processing protocol, or a change in the canopy state around 2010. Analysis of the precipitation records at the tower suggests the first decade was wetter than the second. The average total precipitation for the months of June and July between 2001 and 2010 was about 250 mm, while it was about 200 mm between 2011 and 2020. Since the model does not fully consider surface evaporation following rain events, it could, at least partly, explain the better fit between model and measurement for the drier of the two decades. FIGURE 7. Open in a new tab Diurnal patterns of measured (solid lines) and modeled (dashed lines) latent and sensible heat, as well as the sum of heat terms (latent and sensible heat, soil heat flux for measurements and model, with the addition of stem and leaf heat storage for model). Dotted lines refer to the standard deviation across years. FIGURE 8. Open in a new tab Comparison between measured and modeled sensible heat fluxes at Morgan Monroe between 2000 and 2010 (left), and between 2011 and 2020 (right). The dashed line is the 1:1 relation; the solid line is the linear fit. The net radiation flux densities measured from the top of the canopy compared with the modelled values are shown in Figure 9 . The net radiation is the available energy driving turbulent heat exchanges and must thus be modelled correctly in order to simulate these exchanges accurately. The results shown here suggest a high degree of agreement between model estimates and the measurements, with the EMS and Bartlett sites having significantly more variability in the level of agreement. Reasons for this contrast with results obtained at other sites are not clear. They could be related to measurement issues linked to the radiation sensor being located closer to the canopy top at these sites than at other sites, which could result in more variability due to the smaller sensor footprint. FIGURE 9. Open in a new tab Comparison between net radiation measured from the towers and net radiation modeled. The modeled net radiation is the balance of the PAR, NIR and longwave radiation incoming minus the outgoing, which is equal to the sum of the net radiation of each canopy element (absorbed shortwave and longwave for leaves, wood and soil). 3.2. CO 2 Flux, Heat Storage and Outgoing Longwave Radiation The net ecosystem exchange (NEE) estimates from the model are compared with the measured values in Figure 10 . The agreement is generally good, with the model accounting for about 80%–85% of the variance in NEE measurements. Again, the Hubbard Brook shows a lower agreement, likely due to the data processing not including filtering of fluxes and u*. The saturation in the model estimates at 6 μmol m −2 s −1 is related to the soil respiration being fixed at this value, which is a best estimate in the current absence in CanVeg2 of a component to simulate soil respiration. FIGURE 10. Open in a new tab Comparison of measured and modeled net CO 2 fluxes (Fco2) at all sites used. Data points are hourly or half‐hourly. The dashed line is the 1:1 relation, the solid line is the linear fit. For the Hubbard Brook site data points for which net radiation is less than 50 W/m 2 are shown in blue. Diurnal patterns for measured and modelled soil heat flux densities are shown in Figure 11 , along with modelled stem heat storage, which is a rarely measured quantity at eddy covariance tower sites. The results show a notable difference between measured and modelled soil heat flux ( R 2 around 0.2, see Figures S3 and S4 in Supporting Information ), and a difference in the moment at which the soil heat flux peaks during the day. At the three sites where this measurement was made, soil heat flux peaks in the afternoon, while the model values peak in the morning. It is reasonable to expect this flux density to peak in the afternoon following the accumulation of radiation absorption and a rise in canopy air temperatures. This suggests an issue with this model component, possibly related to the lack of surface evaporation from the soil‐litter layer, and/or the soil thermal properties allowing too much heat transfer too rapidly. However, the measured values appear low in comparison to measurements made near Oak Ridge, TN, by Baldocchi and Vogel ( 1996 ), who reported average soil heat flux peaking around 30 W/m 2 in July using the average from three heat flux plates. Analysis of the raw measurements used in Baldocchi and Vogel ( 1996 ) indicates that soil heat flux peaks around noon (see Figure S5 ), similarly to the model behaviour, suggesting the soil heat flux measurements at Morgan Monroe, Hubbard Brook and Pasoh may not consider the flux of heat occurring above the heat flux plate typically buried at about 15 cm below the surface. It is generally recognised that soil heat flux measurements suffer from large uncertainties due to such practical challenges, and this should be considered when using these measurements to critically evaluate a model. FIGURE 11. Open in a new tab Comparison of diurnal patterns of measured (solid lines) and modeled (dashed lines) soil heat flux, and diurnal pattern of stem heat storage which is not measured and only shown for model. Note that soil heat flux is not measured at the EMS and Bartlett sites. Leaf heat storage included in the model is not shown since the values are very small relative to other fluxes. The stem heat storage in Figure 11 show midday peak flux densities around 60–80 W m −2 at the four temperate sites, and around 90 W m −2 at the tropical site which has more wood structures exposed to solar radiation due to its complex canopy structure. These values are in line with published values by Wilson et al. ( 2000 ), Gu et al. ( 2007 ) and Haverd et al. ( 2007 ). They are, however, significantly higher than the values reported in Oliphant et al. ( 2004 ) which estimated peak flux densities around 10 W m −2 a the Morgan Monroe site, and Ohkubo et al. ( 2008 ) mentioning a study evaluating the peak flux density of biomass heat storage at 15 W m −2 at the Pasoh site. We suggest that these low values attest to the challenge of upscaling a limited number of bole temperature measurements to derive whole canopy heat storage in stems and branches. The modelling of soil heat storage is highly dependent on the amount of radiation absorbed by the soil layer, the turbulent heat exchanges at the soil surface and the conduction of heat through the soil layers. To evaluate the first component, Figure 12 shows the diurnal patterns for the ratio between net radiation at the forest floor and the net radiation at the top of the canopy. Baldocchi and Vogel ( 1996 ) measured this ratio to be around 10% in a temperate broad‐leaved forest located in Tennessee with a similar structure to the temperate sites used here. The values obtained here are close to the 10% ratio and suggest that the modelled net radiation of the soil surface is reasonable. The issue with soil heat flux will be further analysed in the next section below on scalar vertical profiles. FIGURE 12. Open in a new tab Diurnal patterns of modeled ratio between net radiation at the bottom of the canopy and at the top of the canopy for all sites used. Temperature gradients in the canopy are central to a multilayer canopy model like CanVeg2 and are a challenging element to estimate correctly. Reliable measurements of leaf and stem temperatures in different parts of the canopy are difficult to obtain. The emission of longwave radiation by a canopy is a simple measurement which integrates the temperatures of all canopy elements and can be used to evaluate the whole canopy radiative temperature. The comparison between measured and modelled outgoing longwave values is shown in Figure 13 . The model accounted for about 95% of the variance in measurements and suggests the whole canopy temperature is reasonably well simulated. FIGURE 13. Open in a new tab Comparison of measured and modeled outgoing longwave from the top of the canopy at three sites where data was available. At the EMS site, measurements were not made from the EMS tower but from the HDW tower (Hardwood microclimate and walk‐up tower), located nearby in a canopy with very similar species composition and structure as the EMS site. 3.3. Sensitivity Analysis Results of the sensitivity analysis performed at the Morgan Monroe site are shown in Figure 14 . The model was more strongly sensitive to the stem diameter parameter. This is likely because the radiative forcing is calculated on a per surface area basis, and reducing the surface area when calculating the heat storage results in concentrating a higher radiation load on a smaller stem diameter. All other parameters had relatively small effects on the stem heat storage. These results suggest that accurately establishing these parameter values may not be critically important to estimate stem heat storage using a multilayer canopy model, and that uncertainty in the values used for this study did not significantly affect the magnitude of the heat storage estimates obtained. FIGURE 14. Open in a new tab Sensitivity of the stem heat storage to four model parameters at the Morgan Monroe site between 2016 and 2020: Stem diameter profile, stem boundary layer conductance, stem heat capacity and stem thermal conductivity. Solid line is the model estimates using the default model parameter value, the dashed line using half the default value, and the dash‐dotted line double the default value. 3.4. Scalar Vertical Profiles At the EMS site, measurements of CO 2 concentration, relative humidity and air temperature made from different heights above ground were available. The latter two were also available at the Hubbard Brook site. These measurements are valuable to evaluate the model's ability to estimate the sources and sinks of CO 2 , water vapor and heat in the different canopy layers, as well as its ability to disperse the scalars vertically with time. The comparison between the measured and modelled profiles is shown in Figure 15 . The measured profiles generally show higher scalar concentrations near the ground than the model profiles. This is particularly the case for CO 2 concentration at night, early morning and late evening, which could be indicative of the lack of a soil respiration component in CanVeg2. It could also suggest that the dispersion matrix is overestimating the vertical mixing in the lower layers during stable and low wind velocity conditions. The relative humidity is also underestimated by the model in the lower canopy layers, but in this case the underestimation is more pronounced in the middle of the day when the radiative forcing is high. This suggests the lack of enhanced evaporation from the soil‐litter layer following rain events may be at cause, and here also, the dispersion matrix could be overestimating the vertical mixing. The same model limitation in simulating soil surface evaporation may explain the air temperatures in the lower canopy layers being lower in the measurements than in the simulations. The overestimated soil temperatures (Figure 15 bottom right) could also contribute to higher modelled air temperatures compared with measurements. This higher air temperature in lower layers at night may also partly result from positive sensible heat fluxes in those same layers at night mainly caused by the stems being warmer than the air (see the vertical profiles of turbulent heat fluxes shown in Figure S6 in Supporting Information ). FIGURE 15. Open in a new tab Comparison between measured and modeled vertical profiles for relative humidity (%), air temperature (°C) and soil temperature (°C) at the EMS and Hubbard Brook sites for different times of day, and for CO 2 (ppm) for the EMS site. Profiles include all available measurements for the month of July and for each scalar. Solid lines (and dots for relative humidity at Hubbard Brook) refer to the measurements, and dashed line refer to the model. Colors refer to hour (blue: 4 am, green: 8 am, yellow: 2 pm, red: 8 pm and black: Midnight). To assess whether evaporation from the soil surface following rain events may be a factor in the underestimation of relative humidity in the lower canopy layers, we performed an additional analysis at the Harvard Forest site. We looked through the available rainfall data from 2005 and 2020 to find days in the month of July for which there was a daily cumulative rainfall higher than 15 mm with no rainfall above 15 mm in the 3 days prior. We found 19 such events for which we have measurements of relative humidity 2.5 m above ground. We then computed the average residual between the measured and modelled relative humidity at 2.5 m in the 3 days prior to the rain event and 3 days after the event. We found the residual between modelled and measured relative humidity is about 5% higher after the rain event compared with before (median of 4.5%). This suggests the evaporation from wet surfaces is at least partly responsible for the evaporation underestimation. Determining the magnitude of the effect of surface evaporation on relative humidity and latent heat exchanges requires adding this component in the model, which we intend to do in the next model development iteration. Figure 15 shows in the bottom right the comparison between modelled soil temperature profiles and soil temperature measurements made at different depths. The overestimation of soil temperature is in line with the apparent overestimated soil heat flux discussed in the previous section. There is, however, a similar amplitude between the mid‐afternoon and nighttime/early morning soil temperatures at the 5 cm depth, suggesting the soil heat capacity and thermal conductance values used may be adequate, and pointing to turbulent heat exchange and/or radiative forcing as more likely causes of the overestimated soil temperatures through inaccurate soil albedo and/or underestimated latent heat fluxes between the soil surface and the lower canopy air space. Soil optical properties (reflectances in PAR and NIR) values used are average values and may not represent specific site conditions. Attention was given to soil thermal properties in this study, including its vertical variability, and a sensitivity analysis suggests that adjusting soil thermal properties has limited capacity to reduce heat transfer in the soil layers when using an energy balance approach, as the radiative forcing and turbulent heat exchange have the largest influence on the amount of heat forced down the soil layers. The lack of information on the soil surface moisture content may limit the model's ability to correctly simulate surface evaporation, which leads to overestimation of soil temperature, underestimation of relative humidity in the lower canopy layers, and overestimation of air temperature in the lower canopy layers. The underestimation of CO 2 in lower canopy layers could be caused by a combination of the absence of bole respiration in the model, the underestimation of soil respiration, or overestimated turbulent mixing. To explore the latter possibility and assess the sensitivity of the scalar profiles to the dispersion matrix, we did an additional model run using resistances in the dispersion matrix increased by a factor of four. This resulted in a reasonable match between measured and modelled CO 2 profiles; however, the level of correspondence between measured and modelled profiles for other scalars decreased (see Figure S7 in Supporting Information ). Hence, we conclude that the turbulent mixing alone cannot explain discrepancies between measured and modelled scalar profiles, and that bole respiration, and soil evaporation and respiration pulses may contribute. On this basis, we suggest that the measurement of soil‐litter moisture content at flux tower sites has potential to improve modelling of surface evaporation and soil temperature, and ultimately canopy microclimates. 3.5. Areas of Further CanVeg2 Model Improvements Future versions of CanVeg2 would benefit from the integration of precipitation data recorded at flux tower sites to model topsoil moisture content and evaporation from surfaces, including leaves, branches, stems and soil litter. This could substantially improve simulations of soil temperature, canopy latent heat fluxes and eventually soil respiration through the addition of such model component. Currently, soil respiration is a fixed value which represents an average for a forest canopy, and evaporation from the soil surface uses the soil moisture measured at 15 cm depth and ignores the high soil‐litter moisture content following rain events. This likely leads to lower humidity and higher temperature in the lower canopy air layers, higher soil temperatures through underestimated latent heat exchanges and overestimated soil heat fluxes. Given that the model currently estimates the amount of radiation absorbed by stems and branches, as well as their temperatures, a component to simulate bole respiration would likely benefit model accuracy. Edwards and Hanson ( 1996 ) estimated bole respiration to be about 5%–8% of the total yearly carbon flux near the Oak Ridge eddy covariance tower, in an oak forest with a LAI of about 5. The coefficients of determination ( R 2 ) between measured and modelled sensible heat flux densities shown in Figure 6 fall short of the 0.8 threshold which may be considered a typical figure of merit for tests of biophysical models (Baldocchi and Wilson 2001 ). The improvement of evaporation rates from the soil surface could potentially improve the sensible heat flux accuracies. The simulation of sensible heat flux is a challenging quantity to estimate when using an energy balance modelling approach since any inaccuracy in estimates of boundary layer conductances or latent heat exchange at a given model time step can result in important changes in leaf, stem or soil temperatures, eventually leading to “swings” in sensible heat flux density estimates between model time steps as the temperature stabilises. The use of smaller model time steps and Runge–Kutta numerical methods could reduce this effect (Bonan et al. 2026 ). It is also worth noting that errors in air temperature profile will induce errors in sensible heat flux. This is particularly true in lower canopy layers in the afternoon where the difference between mean air temperature modelled versus observed was about 1.5 degrees (see Figure 15 ). 4. Conclusion We evaluated flux density estimates from the CanVeg2 multilayer canopy model against eddy covariance measurements at five broad‐leaf forest sites with contrasting canopy structures and over long‐time series. We also evaluated the model estimates of scalar vertical profiles at one of the sites. The CanVeg2 model builds on the original CanVeg model of Baldocchi and Harley ( 1995 ) and expands its capabilities by integrating 3D radiative transfer simulations, adding stems as canopy elements to estimate stem heat storage and integrating higher order closure principles to derive wind velocity and variance statistics as a function of canopy structure. The CanVeg2 model structure was here described, and the modelled flux density estimates and scalar vertical profiles compared very favourably to measurements. Model estimates of latent heat flux, net ecosystem exchange rates and outgoing longwave radiation accounted for more than 80% of the variance in measurements, while sensible heat flux and soil heat flux had lower coefficients of determination. The scalar vertical profiles highlighted diverging model estimates from measurements in lower canopy layers, suggesting inaccuracies in the fluxes near the soil surface. Does the evaluation presented conclude the CanVeg2 model is accurate, or fit for purpose? We have highlighted areas of strength, some weaknesses and specific areas with potential for improvement. These improvements include measuring or estimating the soil surface moisture content using rainfall data records and using the measured or estimated soil surface moisture content to calculate soil‐litter respiration and improve soil surface evaporation, as well as using the simulated bole temperatures to calculate bole respiration. Vitally, the results show good progress in modelling—from first principles and tuning only four key physiological parameters by site– all canopy fluxes and within canopy microclimate over several years covering contrasting environmental conditions. The study highlighted potential improvement likely to significantly reduce remaining areas of uncertainty to make it fit for the purpose of testing hypotheses of canopy functioning and modelling responses to changes in climate, thus guiding future model development efforts. We conclude the CanVeg2 model nears suitability to mechanistically extrapolate canopy responses into future environmental conditions from underlying principles and predict reductions in productivity following extreme events like concurrent heat and drought if hydraulic safety parameters and soil matric potential are adequately set. Considering these strengths and weaknesses, the model proves a useful and promising tool to address new science questions and improve our understanding of interactions between canopy functions influencing its response to changing climatic conditions and its ability to sustain extreme conditions. The model's strong point is its multilayer nature, which enables a more distinctive temporal description of the state of different canopy parts that can be exposed to very different levels of hydraulic or heat stress, and thus promises to achieve more realistic predictions of whole canopy response to extreme climatic conditions compared with big‐leaf models, which rely on canopy averages to describe non‐linear processes and usually depend more heavily on empirically derived parameters. Author Contributions Hideki Kobayashi: software, writing – review and editing. Dennis Baldocchi: conceptualization, formal analysis, investigation, methodology, software, validation, writing – review and editing. J. William Munger: data curation, resources, writing – review and editing. Gordon B. Bonan: formal analysis, investigation, methodology, software, writing – review and editing. Tilden P. Meyers: conceptualization, software, writing – review and editing. Martin Béland: conceptualization, formal analysis, funding acquisition, investigation, methodology, software, validation, writing – original draft. Funding This work was supported by Natural Sciences and Engineering Research Council of Canada (ALLRP 590324‐23); National Science Foundation (1852977). Conflicts of Interest The authors declare no conflicts of interest. Supporting information Table S1: Model parameters common to all sites. Table S2: Site‐specific model parameters. Figure S1: Model NEE sensitivity to Vcmax at the Morgan Monroe site between 2016 and 2020. Figure S2: Model Bowen ratio sensitivity to iota at the Morgan Monroe site between 2016 and 2020. Figure S3: Comparison of measured and modelled soil heat flux densities at the Morgan Monroe site. Figure S4: Comparison of measured and modelled soil heat flux densities at the Pasoh site. Figure S5: Diurnal pattern of soil heat flux at the walker branch watershed site in July 1995 (top) and 1996 (bottom). Data collection was presented in Baldocchi and Vogel ( 1996 ). Figure S6: Vertical profiles of turbulent heat fluxes at all sites used. Sensible heat includes sources from leaves and wood. Figure S7: Measured and modelled scalar vertical profiles and soil temperature profile at the EMS site for model run using resistances in the dispersion matrix increased by a factor of 4. GCB-32-e70867-s001.docx (1.7MB, docx) Acknowledgements This material is based upon work supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number ALLRP 590324‐23 and the NSF National Center for Atmospheric Research (NCAR), which is a major facility sponsored by the National Science Foundation (NSF) under Cooperative Agreement No. 1852977. M.B. thanks UC Berkeley and NCAR for hosting sabbatical visits in 2023‐2024 during which time much of this research was carried out, and Jeffrey Wood (U. of Missouri) and Ned Patton (NCAR) for fruitful discussions. M.B. thanks Jack Hastings for help with field work at the New Hampshire sites, and for the acquiring the leaf angle measurements at those sites. Also, Andrew Ouimette and Daniel Evans for help with field work at Bartlett and Hubbard Brook, and Eric Kelsey and Mark Green for providing vertical profiles of scalars at Hubbard Brook. The authors acknowledge the Harvard Forest, Morgan Monroe, Bartlett and Hubbard Brooks AmeriFlux sites for their data records. Funding for AmeriFlux data resources was provided by the U.S. Department of Energy's Office of Science. We acknowledge the Pasoh Forest Reserve FLUXNET site for their data records. Data Availability Statement Ecosystem flux data are available at: https://ameriflux.lbl.gov/ for the Harvard Forest, Morgan Monroe, Bartlett and Hubbard Brook sites, and at https://fluxnet.org for the Pasoh site. The data supporting the findings of this study are available at https://doi.org/10.5281/zenodo.17575875 . References Baldocchi, D.

1991. “On Estimating HNO3 Deposition to a Deciduous Frosts With a Lagrangian Random‐Walk Model.” In Precipitation Scavenging and Atmosphere–Surface Exchange, 2, 1081–1093. CRC Press. [ Google Scholar ] Baldocchi, D. , and Bowling D.. 2003. “Modelling the Discrimination of 13 CO 2 Above and Within a Temperate Broad‐Leaved Forest Canopy on Hourly to Seasonal Time Scales.” Plant, Cell & Environment 26, no. 2: 231–244. [ Google Scholar ] Baldocchi, D. , and Meyers T.. 1989. Turbulence Spectra in a Deciduous Forest. American Meteorological Society. [ Google Scholar ] Baldocchi, D. D. , Fuentes J. D., Bowling D. R., Turnipseed A. A., and Monson R. K.. 1999. “Scaling Isoprene Fluxes From Leaves to Canopies: Test Cases Over a Boreal Aspen and a Mixed Species Temperate Forest.” Journal of Applied Meteorology 38, no. 7: 885–898. [ Google Scholar ] Baldocchi, D. D. , and Harley P. C.. 1995. “Scaling Carbon‐Dioxide and Water‐Vapor Exchange From Leaf to Canopy in a Deciduous Forest .2. Model Testing and Application.” Plant, Cell & Environment 18, no. 10: 1157–1173. [ Google Scholar ] Baldocchi, D. D. , Law B. E., and Anthoni P. M.. 2000. “On Measuring and Modeling Energy Fluxes Above the Floor of a Homogeneous and Heterogeneous Conifer Forest.” Agricultural and Forest Meteorology 102, no. 2–3: 187–206. [ Google Scholar ] Baldocchi, D. D. , and Meyers T. P.. 1988. “Turbulence Structure in a Deciduous Forest.” Boundary‐Layer Meteorology 43, no. 4: 345–364. [ Google Scholar ] Baldocchi, D. D. , and Vogel C. A.. 1996. “Energy and CO2 Flux Densities Above and Below a Temperate Broad‐Leaved Forest and a Boreal Pine Forest.” Tree Physiology 16, no. 1–2: 5–16. [ DOI ] [ PubMed ] [ Google Scholar ] Baldocchi, D. D. , and Wilson K. B.. 2001. “Modeling CO 2 and Water Vapor Exchange of a Temperate Broadleaved Forest Across Hourly to Decadal Time Scales.” Ecological Modelling 142, no. 1–2: 155–184. [ Google Scholar ] Baldocchi, D. D. , Ryu Y., Dechant B., et al. 2020. “Outgoing Near Infrared Radiation From Vegetation Scales With Canopy Photosynthesis Across a Spectrum of Function, Structure, Physiological Capacity and Weather.” Journal of Geophysical Research – Biogeosciences 125: e2019JG005534. [ Google Scholar ] Béland, M.

2025. “Mapping Wood Area in Forests From Ground Lidar and Estimating Their Light Interception Using Radiative Transfer Modeling.” Agricultural and Forest Meteorology 375: 110883. [ Google Scholar ] Béland, M. , and Baldocchi D.. 2020. “Is Foliage Clumping an Outcome of Resource Limitations Within Forests?” Agricultural and Forest Meteorology 295: 108185. [ Google Scholar ] Béland, M. , and Baldocchi D. D.. 2021. “Vertical Structure Heterogeneity in Broadleaf Forests: Effects on Light Interception and Canopy Photosynthesis.” Agricultural and Forest Meteorology 307: 108525. [ Google Scholar ] Béland, M. , Baldocchi D. D., Widlowski J.‐L., Fournier R. A., and Verstraete M. M.. 2014. “On Seeing the Wood From the Leaves and the Role of Voxel Size in Determining Leaf Area Distribution of Forests With Terrestrial LiDAR.” Agricultural and Forest Meteorology 184: 82–97. [ Google Scholar ] Béland, M. , Bonan G. B., Kobayashi H., and Baldocchi D.. in review. “Modifications in the Norman (1979) Canopy Radiative Transfer Model to Consider Multiple Scattering in Clumped Layers and Radiation Absorption by Stems.” Journal of Advances in Modeling Earth Systems. [ Google Scholar ] Béland, M. , and Kobayashi H.. 2021. “Mapping Forest Leaf Area Density From Multiview Terrestrial Lidar.” Methods in Ecology and Evolution 12, no. 4: 619–633. [ Google Scholar ] Béland, M. , and Kobayashi H.. 2024. “Drivers of Deciduous Forest Near‐Infrared Reflectance: A 3D Radiative Transfer Modeling Exercise Based on Ground Lidar.” Remote Sensing of Environment 302: 113951. [ Google Scholar ] Béland, M. , Widlowski J.‐L., Fournier R., Côté J.‐F., and Verstraete M. M.. 2011. “Estimating Leaf Area Distribution in Savanna Trees From Terrestrial LiDAR Measurements.” Agricultural and Forest Meteorology 151, no. 9: 1252–1266. [ Google Scholar ] Beyschlag, W. , and Ryel R. J.. 2007. Canopy Photosynthesis Modeling, Functional Plant Ecology. CRC Press. [ Google Scholar ] Bonan, G.

2019. Climate Change and Terrestrial Ecosystem Modeling. Cambridge University Press. [ Google Scholar ] Bonan, G. B. , Burns S. P., and Patton E. G.. 2026. “Beyond Surface Fluxes: Observational and Computational Needs of Multilayer Canopy Models—A Walnut Orchard Test Case.” Agricultural and Forest Meteorology 378: 110960. [ Google Scholar ] Bonan, G. B. , Patton E. G., Finnigan J. J., Baldocchi D. D., and Harman I. N.. 2021. “Moving Beyond the Incorrect but Useful Paradigm: Reevaluating Big‐Leaf and Multilayer Plant Canopies to Model Biosphere‐Atmosphere Fluxes—A Review.” Agricultural and Forest Meteorology 306: 108435. [ Google Scholar ] Bonan, G. B. , Patton E. G., Harman I. N., et al. 2018. “Modeling Canopy‐Induced Turbulence in the Earth System: A Unified Parameterization of Turbulent Exchange Within Plant Canopies and the Roughness Sublayer (CLM‐ml v0).” Geoscientific Model Development 11, no. 4: 1467–1496. [ Google Scholar ] Bonan, G. B. , Williams M., Fisher R. A., and Oleson K. W.. 2014. “Modeling Stomatal Conductance in the Earth System: Linking Leaf Water‐Use Efficiency and Water Transport Along the Soil–Plant–Atmosphere Continuum.” Geoscientific Model Development 7, no. 5: 2193–2222. [ Google Scholar ] Campbell, G. S.

1974. “A Simple Method for Determining Unsaturated Conductivity From Moisture Retention Data.” Soil Science 117, no. 6: 311–314. [ Google Scholar ] Choi, M. , Jacobs J. M., and Kustas W. P.. 2008. “Assessment of Clear and Cloudy Sky Parameterizations for Daily Downwelling Longwave Radiation Over Different Land Surfaces in Florida, USA.” Geophysical Research Letters 35, no. 20: 2008GL035731. [ Google Scholar ] Christoffersen, B. O. , Gloor M., Fauset S., et al. 2016. “Linking Hydraulic Traits to Tropical Forest Function in a Size‐Structured and Trait‐Driven Model (TFS v.1‐Hydro).” Geoscientific Model Development 9, no. 11: 4227–4255. [ Google Scholar ] Collatz, G. , Berry J., Farquhar G., and Pierce J.. 1990. “The Relationship Between the Rubisco Reaction Mechanism and Models of Photosynthesis.” Plant, Cell & Environment 13, no. 3: 219–225. [ Google Scholar ] Cowan, I.

1965. “Transport of Water in the Soil‐Plant‐Atmosphere System.” Journal of Applied Ecology 2: 221–239. [ Google Scholar ] Cowan, I.

1978. “Stomatal Behaviour and Environment.” Advances in Botanical Research 4: 117–228. [ Google Scholar ] Daamen, C. C. , and Simmonds L. P.. 1996. “Measurement of Evaporation From Bare Soil and Its Estimation Using Surface Resistance.” Water Resources Research 32, no. 5: 1393–1402. [ Google Scholar ] dePury, D. G. G. , and Farquhar G. D.. 1997. “Simple Scaling of Photosynthesis From Leaves to Canopies Without the Errors of Big‐Leaf Models.” Plant, Cell & Environment 20, no. 5: 537–557. [ Google Scholar ] Edwards, N. T. , and Hanson P. J.. 1996. “Stem Respiration in a Closed‐Canopy Upland Oak Forest.” Tree Physiology 16, no. 4: 433–439. [ DOI ] [ PubMed ] [ Google Scholar ] Farquhar, G. D. , von Caemmerer S. V., and Berry J. A.. 1980. “A Biochemical Model of Photosynthetic CO 2 Assimilation in Leaves of C3 Species.” Planta 149, no. 1: 78–90. [ DOI ] [ PubMed ] [ Google Scholar ] Goudriaan, J.

1977. Crop Micrometeorology: A Simulation Study. Wageningen University and Research. [ Google Scholar ] Grace, J. , and Wilson J.. 1976. “The Boundary Layer Over a Populus Leaf.” Journal of Experimental Botany 27, no. 2: 231–241. [ Google Scholar ] Gu, L. , Meyers T., Pallardy S. G., et al. 2007. “Influences of Biomass Heat and Biochemical Energy Storages on the Land Surface Fluxes and Radiative Temperature.” Journal of Geophysical Research 112, no. D2: D02107. [ Google Scholar ] Harley, P. C. , and Baldocchi D. D.. 1995. “Scaling Carbon Dioxide and Water Vapour Exchange From Leaf to Canopy in a Deciduous Forest. I. Leaf Model Parametrization.” Plant, Cell and Environment 18, no. 10: 1146–1156. [ Google Scholar ] Hastings, J.

2025. “Personnal Communication.” Haverd, V. , Cuntz M., Leuning R., and Keith H.. 2007. “Air and Biomass Heat Storage Fluxes in a Forest Canopy: Calculation Within a Soil Vegetation Atmosphere Transfer Model.” Agricultural and Forest Meteorology 147, no. 3–4: 125–139. [ Google Scholar ] Helliker, B. R. , Song X., Goulden M. L., et al. 2018. “Assessing the Interplay Between Canopy Energy Balance and Photosynthesis With Cellulose Delta(18)O: Large‐Scale Patterns and Independent Ground‐Truthing.” Oecologia 187, no. 4: 995–1007. [ DOI ] [ PubMed ] [ Google Scholar ] Jarvis, P. , Miranda H., and Muetzelfeldt R.. 1985. “Modelling Canopy Exchanges of Water Vapor and Carbon Dioxide in Coniferous Forest Plantations.” In the Forest‐Atmosphere Interaction: Proceedings of the Forest Environmental Measurements Conference Held at Oak Ridge, Tennessee, October 23–28, 1983, 521–542. Springer. [ Google Scholar ] Kelsey, E. , and Green M.. 2025. “AmeriFlux BASE US‐HBK Hubbard Brook Experimental Forest, Ver. 3‐5, AmeriFlux AMP, (Dataset).” 10.17190/AMF/1634881. [ DOI ] Kira, T. , Manokaran N., and Appanah S.. 2013. NPP Tropical Forest: Pasoh, Malaysia, 1971–1973, R1. ORNL DAAC. [ Google Scholar ] Knohl, A. , and Baldocchi D. D.. 2008. “Effects of Diffuse Radiation on Canopy Gas Exchange Processes in a Forest Ecosystem.” Journal of Geophysical Research: Biogeosciences 113, no. G2: G02023. [ Google Scholar ] Kobayashi, H. , Baldocchi D. D., Ryu Y., et al. 2012. “Modeling Energy and Carbon Fluxes in a Heterogeneous Oak Woodland: A Three‐Dimensional Approach.” Agricultural and Forest Meteorology 152, no. 2: 83–100. [ Google Scholar ] Kobayashi, H. , and Iwabuchi H.. 2008. “A Coupled 1‐D Atmosphere and 3‐D Canopy Radiative Transfer Model for Canopy Reflectance, Light Environment, and Photosynthesis Simulation in a Heterogeneous Landscape.” Remote Sensing of Environment 112, no. 1: 173–185. [ Google Scholar ] Kosugi, Y. , and Takanashi S.. 2019. “Personal Communication.” Landsberg, J. , and Powell D.. 1973. “Surface Exchange Characteristics of Leaves Subject to Mutual Interference.” Agricultural Meteorology 12: 169–184. [ Google Scholar ] Lawrence, D. M. , Fisher R. A., Koven C. D., et al. 2019. “The Community Land Model Version 5: Description of New Features, Benchmarking, and Impact of Forcing Uncertainty.” Journal of Advances in Modeling Earth Systems 11, no. 12: 4245–4287. [ Google Scholar ] Majasalmi, T. , and Bright R.. 2019. “Evaluation of Leaf‐Level Optical Properties Employed in Land Surface Models—Example With CLM 5.0.” Geoscientific Model Development Discussions 12, no. 9: 3923–3938. [ Google Scholar ] Matthes, J. H. , Munger J. W., and Wofsy S.. 2025. “Biomass Inventories at Harvard Forest EMS Tower Since 1993 ver 44.” Environmental Data Initiative. 10.6073/pasta/c6c57614ba20dd50cee401ec700d3913. [ DOI ] Medlyn, B.

2004. “A maestro retrospective.” In Forests at the Land–Atmosphere Interface, 105–121. CAB International. [ Google Scholar ] Medlyn, B. E. , Robinson A. P., Clement R., and McMurtrie R. E.. 2005. “On the Validation of Models of Forest CO 2 Exchange Using Eddy Covariance Data: Some Perils and Pitfalls.” Tree Physiology 25, no. 7: 839–857. [ DOI ] [ PubMed ] [ Google Scholar ] Meyers, T. , and Tha Paw U K.. 1986. “Testing of a Higher‐Order Closure Model for Modeling Airflow Within and Above Plant Canopies.” Boundary‐Layer Meteorology 37, no. 3: 297–311. [ Google Scholar ] Meyers, T. P. , and Tha Paw U K.. 1987. “Modelling the Plant Canopy Micrometeorology With Higher‐Order Closure Principles.” Agricultural and Forest Meteorology 41, no. 1–2: 143–163. [ Google Scholar ] Monteith, J. , and Unsworth M.. 2013. Principles of Environmental Physics: Plants, Animals, and the Atmosphere. Academic Press. [ Google Scholar ] Munger, J. W.

2022. AmeriFlux FLUXNET‐1F US‐Ha1 Harvard Forest EMS Tower (HFR1), Ver. 3‐5, AmeriFlux AMP, (Dataset). Munger, W. , and Wofsy S.. 2024. Canopy‐Atmosphere Exchange of Carbon, Water and Energy at Harvard Forest EMS Tower since 1991 ver 36. Environmental Data Initiative. 10.6073/pasta/56c6fe02a07e8a8aaff44a43a9d9a6a5. [ DOI ] Norman, J.

1979. “Modeling the Complete Crop Canopy.” In Modification of the Aerial Environment of Crops, edited by EM11 , 249–277. American Society Agricultural Engineers. [ Google Scholar ] Novick, K. , and Phillips R.. 2022. AmeriFlux FLUXNET‐1F US‐MMS Morgan Monroe State Forest, Ver. 3‐5, AmeriFlux AMP, (Dataset). Ohkubo, S. , Kosugi Y., Takanashi S., Matsuo N., Tani M., and Nik A. R.. 2008. “Vertical Profiles and Storage Fluxes of CO 2 , Heat and Water in a Tropical Rainforest at Pasoh, Peninsular Malaysia.” Tellus Series B: Chemical and Physical Meteorology 60, no. 4: 569–582. [ Google Scholar ] Oliphant, A. J. , Grimmond C. S. B., Zutter H. N., et al. 2004. “Heat Storage and Energy Balance Fluxes for a Temperate Deciduous Forest.” Agricultural and Forest Meteorology 126, no. 3–4: 185–201. [ Google Scholar ] Oliphant, A. J. , and Stoy P. C.. 2018. “An Evaluation of Semiempirical Models for Partitioning Photosynthetically Active Radiation Into Diffuse and Direct Beam Components.” Journal of Geophysical Research: Biogeosciences 123, no. 3: 889–901. [ Google Scholar ] Ouimette, A. P. , Ollinger S. V., Richardson A. D., et al. 2018. “Carbon Fluxes and Interannual Drivers in a Temperate Forest Ecosystem Assessed Through Comparison of Top‐Down and Bottom‐Up Approaches.” Agricultural and Forest Meteorology 256‐257: 420–430. [ Google Scholar ] Pastorello, G. , Trotta C., Canfora E., et al. 2020. “The FLUXNET2015 Dataset and the ONEFlux Processing Pipeline for Eddy Covariance Data.” Sci Data 7, no. 1: 225. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Pyles, R. D. , Weare B. C., Tha Paw U K., and Gustafson W.. 2003. “Coupling Between the University of California, Davis, Advanced Canopy–Atmosphere–Soil Algorithm (ACASA) and MM5: Preliminary Results for July 1998 for Western North America.” Journal of Applied Meteorology 42, no. 5: 557–569. [ Google Scholar ] Rastetter, E. B. , King A. W., Cosby B. J., Hornberger G. M., O'Neill R. V., and Hobbie J. E.. 1992. “Aggregating Fine‐Scale Ecological Knowledge to Model Coarser‐Scale Attributes of Ecosystems.” Ecological Applications 2, no. 1: 55–70. [ DOI ] [ PubMed ] [ Google Scholar ] Raupach, M.

1988. Canopy Transport Processes, Flow and Transport in the Natural Environment: Advances and Applications. Springer. [ Google Scholar ] Raupach, M. R.

1989. “A Practical Lagrangian Method for Relating Scalar Concentrations to Source Distributions in Vegetation Canopies.” Quarterly Journal of the Royal Meteorological Society 115, no. 487: 609–632. [ Google Scholar ] Raupach, M. R. , and Finnigan J. J.. 1988. “Single‐Layer Models of Evaporation From Plant Canopies Are Incorrect but Useful, Whereas Multilayer Models Are Correct but Useless—Discuss.” Australian Journal of Plant Physiology 15, no. 6: 705–716. [ Google Scholar ] Ryu, Y. , Sonnentag O., Nilson T., et al. 2010. “How to Quantify Tree Leaf Area Index in an Open Savanna Ecosystem: A Multi‐Instrument and Multi‐Model Approach.” Agricultural and Forest Meteorology 150, no. 1: 63–76. [ Google Scholar ] Schuepp, P.

1993. “Tansley Review No. 59. Leaf Boundary Layers.” New Phytologist 125: 477–507. [ DOI ] [ PubMed ] [ Google Scholar ] Sellers, P. , Randall D. A., Collatz G. J., et al. 1996. “A Revised Land Surface Parameterization (SiB2) for Atmospheric GCMs. Part I: Model Formulation.” Journal of Climate 9, no. 4: 676–705. [ Google Scholar ] Swenson, S. C. , Burns S. P., and Lawrence D. M.. 2019. “The Impact of Biomass Heat Storage on the Canopy Energy Balance and Atmospheric Stability in the Community Land Model.” Journal of Advances in Modeling Earth Systems 11, no. 1: 83–98. [ Google Scholar ] Waggoner, P. E. , and Reifsnyder W. E.. 1968. “Simulation of the Temperature, Humidity and Evaporation Profiles in a Leaf Canopy.” Journal of Applied Meteorology and Climatology 7, no. 3: 400–409. [ Google Scholar ] Wang, Y. , and Jarvis P.. 1990. “Description and Validation of an Array Model—MAESTRO.” Agricultural and Forest Meteorology 51, no. 3: 257–280. [ Google Scholar ] Wang, Y.‐P. , and Leuning R.. 1998. “A Two‐Leaf Model for Canopy Conductance, Photosynthesis and Partitioning of Available Energy I: Model Description and Comparison With a Multi‐Layered Model.” Agricultural and Forest Meteorology 91, no. 1–2: 89–111. [ Google Scholar ] Wayson, C. A. , Randolph J. C., Hanson P. J., Grimmond C. S. B., and Schmid H. P.. 2006. “Comparison of Soil Respiration Methods in a Mid‐Latitude Deciduous Forest.” Biogeochemistry 80, no. 2: 173–189. [ Google Scholar ] Weiss, A. , and Norman J.. 1985. “Partitioning Solar Radiation Into Direct and Diffuse, Visible and Near‐Infrared Components.” Agricultural and Forest Meteorology 34, no. 2–3: 205–213. [ Google Scholar ] Williams, M. , Rastetter E. B., Fernandes D. N., et al. 1996. “Modelling the Soil‐Plant‐Atmosphere Continuum in a Quercus‐Acer Stand at Harvard Forest: The Regulation of Stomatal Conductance by Light, Nitrogen and Soil/Plant Hydraulic Properties.” Plant, Cell & Environment 19, no. 8: 911–927. [ Google Scholar ] Wilson, K. , Hanson P., and Baldocchi D.. 2000. “Factors Controlling Evaporation and Energy Partitioning Beneath a Deciduous Forest Over an Annual Cycle.pdf.” Wilson, T. , Norman J., Bland W., and Kucharik C.. 2003. “Evaluation of the Importance of Lagrangian Canopy Turbulence Formulations in a Soil–Plant–Atmosphere Model.” Agricultural and Forest Meteorology 115, no. 1–2: 51–69. [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Table S1: Model parameters common to all sites. Table S2: Site‐specific model parameters. Figure S1: Model NEE sensitivity to Vcmax at the Morgan Monroe site between 2016 and 2020. Figure S2: Model Bowen ratio sensitivity to iota at the Morgan Monroe site between 2016 and 2020. Figure S3: Comparison of measured and modelled soil heat flux densities at the Morgan Monroe site. Figure S4: Comparison of measured and modelled soil heat flux densities at the Pasoh site. Figure S5: Diurnal pattern of soil heat flux at the walker branch watershed site in July 1995 (top) and 1996 (bottom). Data collection was presented in Baldocchi and Vogel ( 1996 ). Figure S6: Vertical profiles of turbulent heat fluxes at all sites used. Sensible heat includes sources from leaves and wood. Figure S7: Measured and modelled scalar vertical profiles and soil temperature profile at the EMS site for model run using resistances in the dispersion matrix increased by a factor of 4. GCB-32-e70867-s001.docx (1.7MB, docx) Data Availability Statement Ecosystem flux data are available at: https://ameriflux.lbl.gov/ for the Harvard Forest, Morgan Monroe, Bartlett and Hubbard Brook sites, and at https://fluxnet.org for the Pasoh site. The data supporting the findings of this study are available at https://doi.org/10.5281/zenodo.17575875 . Articles from Global Change Biology are provided here courtesy of Wiley ACTIONS View on publisher site PDF (9.4 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top