Modeling photosynthesis of discontinuous plant canopies by linking the Geometric Optical Radiative Transfer model with biochemical processes

Modeling vegetation photosynthesis is essential for understanding carbon exchanges between terrestrial ecosystems and the atmosphere. The radiative transfer process within plant canopies is one of the key drivers that regulate canopy photosynthesis. Most vegetation cover consists of discrete plant crowns, of which the physical observation departs from the underlying assumption of a homogenous and uniform medium in classic radiative transfer theory. Here we advance the Geometric Optical Radiative Transfer (GORT) model to simulate photosynthesis activities for discontinuous plant canopies. We separate radiation absorption into two components that are absorbed by sunlit and shaded leaves, and derive analytical solutions by integrating over the canopy layer. To model leaf-level and canopylevel photosynthesis, leaf light absorption is then linked to the biochemical process of gas diffusion through leaf stomata. The canopy gap probability derived from GORT differs from classic radiative transfer theory, especially when the leaf area index is high, due to leaf clumping effects. Tree characteristics such as tree density, crown shape, and canopy length affect leaf clumping and regulate radiation interception. Modeled gross primary production (GPP) for two deciduous forest stands could explain more than 80 % of the variance of flux tower measurements at both near hourly and daily timescales. We demonstrate that ambient CO2 concentrations influence daytime vegetation photosynthesis, which needs to be considered in biogeochemical models. The proposed model is complementary to classic radiative transfer theory and shows promise in modeling the radiative transfer process and photosynthetic activities over discontinuous forest canopies.


Introduction
Terrestrial plants assimilate atmospheric carbon dioxide through photosynthesis (Keenan et al., 2013;Myneni et al., 1997).The climate system, in turn, affects vegetation development and photosynthetic activities (Broich et al., 2014;Xia et al., 2014;Yi et al., 2010).Photosynthesis, accompanied by exchanges of heat, water vapor, and trace gases within the planetary boundary layer, modifies microclimates and local environments and determines ecosystem functions and services (Peng et al., 2014;Xu et al., 2013).The complex biosphere/atmosphere feedbacks are dynamic and interactive (Bonan, 2008;Heimann and Reichstein, 2008), such that robust numerical models that simulate vegetation photosynthesis are required in terrestrial ecosystem models to understand the global carbon cycle (Cramer et al., 2001;Kucharik et al., 2006).
Vegetation photosynthesis activity is regulated by environmental factors, and the light environment within plant canopies is one of the key drivers (Law et al., 2002;Pearcy and Sims, 1994).Biophysical models such as production efficiency models assume linear relationships between absorbed photosynthetically active radiation (APAR) and vegetation primary production (Field et al., 1995;Monteith, 1977;Potter et al., 1993;Prince and Goward, 1995;Running et Q.Xin et al.: Modeling photosynthesis of discontinuous plant canopies al., 2000).Because vegetation photosynthesis harvests solar radiation by green chlorophyll, recent studies have attempted to quantify the fractions of APAR that are absorbed by green chlorophyll (Zhang et al., 2005(Zhang et al., , 2014)).Physiologically, plants assimilate carbon dioxide via the biochemical diffusion processes through stomata, numerous small pores on the leaf surfaces (Collatz et al., 1991;Farquhar and Sharkey, 1982).Stomata can open and close in response to microenvironments, thereby regulating plant carbon uptake (Bonan, 2002).Field physiological studies have accumulated detailed information on the behavior of stomata under certain environmental conditions (Schulze et al., 1994), in which sunlight irradiance plays a vital role (Ball et al., 1987).In this domain, linking the physical process of radiative transfer within plant canopies with the biochemical process of gas diffusion through leaf stomata is essential for accurate representation of vegetation photosynthesis.
Radiative transfer within a plant canopy is determined by many factors such as the partition of incoming solar radiation, solar illumination geometry, terrain slope and aspects, canopy structure, leaf angle distribution, and leaf and substrate spectral properties (Baldocchi et al., 1985;Fan et al., 2014;Schaaf et al., 1994).Classic radiative transfer theory assumes that plant leaves are randomly distributed in three-dimensional space within a homogeneous canopy layer (Goudriaan, 1977;Myneni et al., 1990).The canopy radiative transfer process can be simply characterized by leaf area index (LAI) and leaf angle distribution (LAD).Three-dimensional, multi-layer, and two-leaf radiative transfer models have been developed to simulate leaf absorption of solar irradiance and canopy photosynthesis (Myneni, 1991;Pury and Farquhar, 1997;Ryu et al., 2011;Sellers, 1985).Although classic radiative transfer theory holds well for dense vegetation canopies, most vegetation canopies, especially arboreal canopies, consist of discrete crowns in reality (Yuan et al., 2013).Leaves are clumped within individual crowns, such that more sunlight penetrates to understory layers and the ground surfaces (He et al., 2012;Ni-Meister et al., 2010).Tree crowns also cast shadows on one another and on the background, resulting in self-shadowing effects as described by the geometric-optical theory (Li and Strahler, 1992).Given natural differences in the radiative transfer process between homogenous and discontinuous plant canopies, it is important to understand and account for the influence of crown shape and tree structure on canopy radiation absorption and vegetation photosynthesis.
To address the radiative transfer process in discontinuous canopies, the Geometric Optical Radiative Transfer (GORT) model conceptually combines geometric optical principles for canopy structure and radiative transfer theory for volumetric scattering within canopy crowns (Li et al., 1995).The geometric optical method is used to characterize the process by which sunlight passes directly to the ground surface without reaching any canopy crowns.The radiative transfer principle is applied to model the probability of light pen-etration as it travels through crowns in the canopy.GORT has been used to model the physical aspects of discontinuous plant canopies such as gap fraction, radiation transmission, and bi-directional reflectance (Ni et al., 1997(Ni et al., , 1999;;Xin et al., 2012), and has been validated under a variety of environmental conditions (Liu et al., 2008).Recent efforts have been made to develop and evaluate a simplified GORT model for the use in coupled global dynamic terrestrial ecosystem models (Ni-Meister et al., 2010;Yang et al., 2010).Despite these successful applications, the current version of the GORT model does not have analytical solutions for radiation absorption by sunlit and shaded leaves, though previous studies have tried to solve the process of multiple scattering between canopy and background in an iterative manner (Song et al., 2009).However, sunlit and shaded leaves must be treated separately in photosynthesis modeling because flux densities of photosynthetically active radiation (PAR) incident on leaf surfaces are different (He et al., 2013).It is also necessary to integrate vertically over the canopy to derive mean PAR absorbed by sunlit and shaded leaves because of the non-linear light attenuation within the canopy and the non-linear dependence of leaf stomatal conductance on light absorption (Campbell and Norman, 1998).
The objectives of this study are to (1) advance the GORT model by providing analytical solutions to the radiation absorption of sunlit and shaded leaves and (2) link the radiative transfer process to biochemical processes to simulate leaf and canopy photosynthesis.We first describe the principles of our model and then perform model validation with eddy covariance data from two flux towers situated in the New England region of the United States.

Brief description of canopy gap probability modeled using GORT
Gap probability, the probability of photons reaching a given canopy depth without being intercepted by canopy elements, is key to characterizing the radiation distribution within plant canopies.A detailed description for modeling the gap probability with GORT is described in previous studies (Li et al., 1995;Ni et al., 1999); we summarize it briefly here because the concept of gap probability is necessary for understanding our subsequent work.For homogeneous canopies, Beer's law describes the gap probability of sunlight penetration.For discontinuous plant canopies, leaves are clumped within individual canopy crowns, forming an uneven distribution of gap probabilities for beam radiation.GORT models tree crowns as a collection of ellipsoids (Fig. 1), of which the centers are randomly distributed between the upper and lower boundaries of the canopy layer (h 1 and h 2 ).Each ellipsoid, or each canopy crown, is characterized by one-half of the vertical crown length (b) and a horizontal crown radius (R).The total gap probability is modeled separately as the proportion of sunlight passing through the canopy layer without reaching any crown (hereafter referred to as between-crown gaps) and the proportion of sunlight passing through crowns without being intercepted by canopy leaves (hereafter referred to as withincrown gaps), such that where P gap (h, θ i ) is the gap probability for beam radiation at height h given an illumination zenith angle θ i , P gap (n = 0|h, θ i ) is the between-crown gap, and P gap (n > 0|h, θ i ) is the within-crown gap.
The between-crown gap is modeled based on Boolean theory as an exponential function of crown numbers within a geometric volume that contains no crown centers: where λ v is the tree density, and V is the beam projected cylinder volume with a radius R starting from the canopy top and extending to height h.Assuming that leaves are randomly distributed within each individual crown, the within-crown gap is modeled based on Beer's law as light penetration along the traveling path length, such that where τ (θ i , α) = k b (θ i , α) • FAVD, FAVD is the foliage area volume density within a single crown, and k b (θ i , α) is the extinction coefficient for beam radiation given a specific solar illumination angle θ i and leaf distribution angle α.For a spherical leaf angle distribution, k b = 0.5 cos(θ i ) .P (s|h, θ i ) is the probability distribution function associated with withincrown path length s.
The probability distribution of within-crown paths length can be solved in a convolutional manner: where P (s|n, z, h, θ i ) is the probability distribution of within-crown path length given that a solar ray enters the crown at height h and angle θ i , and P (n|z, h, θ i ) is the probability distribution of the numbers of crowns intercepted by the solar ray incident at angle θ i , entering crowns at height z, and then traveling to height h.Diffuse radiation (i.e., the hemispherically isotropic radiation) can be treated as beam radiation from all directions in the upper hemisphere.The "openness" of discontinuous plant canopies to diffuse radiation on a horizontal plane is defined as where K open (n = 0|h) and K open (n > 0|h) are betweencrown and within-crown openness factors, respectively.θ i is the solar illumination angle, and φ is the azimuth angle.

Sunlit and shaded leaf area index
The gap probability describes the probability of beam radiation being intercepted by plant leaves and hence determines the proportion of leaf areas that are sunlit.For a very thin layer, the reduction of total gap probability is due to leaf interception, the process of which still follows Beer's law: where k b is the canopy extinction coefficient for beam irradiance, δLAI (h) is the leaf area index within a thin layer δh at height h, and P gap (h, θ i ) is the gap probability modeled using GORT.

Q. Xin et al.: Modeling photosynthesis of discontinuous plant canopies
In the limit, as δh becomes infinitely small, we have where P gap (h, θ i ) is the first derivative of gap probability P gap (h, θ i ) with respect to height h.Combining Eqs. ( 5), (9), and (10), we obtain For diffuse radiation, it can be derived in a similar manner: where k d is the extinction coefficient for diffuse irradiance, and K open (h) is the first derivative of the openness factor K open (h) with respect to height h.
The sunlit LAI at height h is the product of the probability of beam sunlight penetration to height h and the probability of sunlight being intercepted by the thin layer and divided by the ratio of leaf area projected on a horizontal surface (Campbell and Norman, 1998), such that where δLAI Sun (h, θ i ) is the sunlit leaf area index within a thin layer δh at height h.Substituting Eqs. ( 9) and (6) into Eq.( 8), we obtain Sunlit LAI for the entire canopy at zenith angle θ is then obtained by integrating from the canopy top to canopy bottom, such that where P gap (h = z 2 |θ i ) and P gap (h = z 1 |θ i ) are the gap probabilities at the canopy top z 2 and canopy bottom z 1 , respectively, whereas the gap probability at the canopy top is 1.
It is worth noting that our calculation of sunlit leaf area for discontinuous canopies is analogous to that for homogeneous canopies, which is given as where LAI * Sun (θ i ) is the sunlit leaf area for homogeneous canopies.
The shaded LAI is simply the remainder of the canopy LAI: (17)

Analytical solutions for the scattering parameters of discontinuous canopies
Canopy scattering parameters such as directionalhemispherical reflectance and hemispherical-hemispherical reflectance (or black-sky albedo and white-sky albedo, respectively) can be obtained by resolving the radiative transfer process or can be approximated using simple analytical solutions.For semi-infinite horizontally homogeneous media, Hapke's solutions of the proportion of unintercepted direct beam (t 0 (h, θ , and directional-hemispherical transmittance (T ∞ df ) are given as (Hapke, 1981) where σ is the single scattering albedo, τ = k(θ i ) L e H is the projected foliage area volume density for the plant canopy, L e is the effective leaf area index, H is the depth of the canopy, θ i is the solar illumination angle, µ i = cos(θ i ) and γ = √ 1 − σ .Starting with surface energy balances, Ni (1998) derived the scattering parameters for a horizontally homogeneous canopy layer with finite thickness as where t ff (h), ρ ff (h), t df (h, θ i ), and ρ df (h, θ i ) are hemispherical-hemispherical transmittance, hemisphericalhemispherical reflectance, directional-hemispherical transmittance, and directional-hemispherical reflectance, respectively.
The scattering parameters for a discontinuous canopy can then be approximated as combinations of a homogeneous where t ff (h), ρ ff (h), t df (h, θ i ), and ρ df (h, θ i ) are hemispherical-hemispherical transmittance, hemisphericalhemispherical reflectance, directional-hemispherical transmittance, and directional-hemispherical reflectance, respectively.Note that our equations here are slightly different from those used by Ni et al. (1999) because between-crown gaps, within which light attenuation obeys Beer's law, are considered in the homogeneous vegetation layer.
The analytical approximation of the canopy reflectance for beam and diffuse radiation is the sum of three factors in radiative transfer: the incoming irradiance scattered by the canopy elements, the first-order scattered radiation from soil background, and the irradiance scattered back and forth between the canopy layer and background surface (Ni et al., 1999).Taking beam radiation as an example and assuming that the background surface is Lambertian, the incoming irradiance scattered by the canopy elements is ρ df , the first-order scattered radiance from soil background is t df ρ s t ff , and the multiple scattering between the canopy elements and soil background is The canopy reflectance for beam irradiance can then be written as The canopy reflectance for diffuse irradiance can be obtained similarly as

Mean photosynthetically active radiation absorbed by sunlit and shaded leaves
Let I 0 be the flux density of incoming solar radiation on a horizontal plane at the top of the canopy and f b be the fraction of incident beam radiation, the unintercepted beam and diffuse fluxes are then where ρ cb and ρ cd are canopy reflectance for beam and diffuse irradiance, respectively; I b and I d are the unintercepted beam and diffuse fluxes, respectively; and k b and k d are canopy extinction coefficients for beam and diffuse irradiance, respectively.The downward beam flux I b is derived based on the assumption of black leaves, meaning that leaves absorb incident irradiance completely and do not transmit radiation (Bonan, 2002).To account for the effects of leaf scattering, the total beam I bt (i.e., unintercepted beam and down scattered beam) and total diffuse I dt (i.e., unintercepted diffuse and down scattered diffuse) irradiance can be modeled by introducing a factor of √ 1 − σ to extinction coefficients similar to the two-stream radiative transfer model (Sellers, 1985).As single scattering albedo increases, the effective extinction coefficient becomes smaller and more sunlight is allowed to transmit through the canopy.That is, where σ is the single scattering albedo of leaves.σ = ρ l + t l , where ρ l and t l are leaf reflectance and transmittance, respectively.
The total irradiance absorbed by the entire canopy per unit ground area consists of leaf absorption for both beam and diffuse irradiance: Substituting Eqs. ( 11), ( 12), (35), and (36) into Eq.( 37), we have Irradiance absorbed by sunlit leaves per unit ground area is obtained as the sum of direct beam, downward scattered beam, and diffuse components:

Q. Xin et al.: Modeling photosynthesis of discontinuous plant canopies
Combining Eqs. ( 33), ( 35), (36), and (40), we have Note that σ is used instead of ρ cd for the beam irradiance of sunlit leaves because sunlit leaves scatter direct beam sunlight only once.
The irradiance absorbed by shaded leaves per unit ground area is simply the difference between the total irradiance absorbed by the canopy and the irradiance absorbed by sunlit leaves: The mean absorbed irradiance for sunlit and shaded canopy per leaf hemi-surface area is then

Modeling leaf photosynthesis and scaling up to canopy photosynthesis
The biochemical process of carbon dioxide assimilation by leaves can be considered as a gas diffusion process through stomata.According to Fick's law, the process is described as where A is the CO 2 assimilation rate, g c is the stomatal conductance, and C a and C i are ambient and intercellular CO 2 concentrations, respectively.Field studies have firmly established the relationship between leaf stomatal conductance and environmental conditions.Jarvis and McNaughton (1986) successfully synthesize the response functions in a multiple-constraint model: where g cmax is the maximum leaf stomatal conductance when environmental factors do not limit carbon uptake and f (x i ) are scalars that account for the influences of various environmental stresses on leaf stomatal conductance.Different formulas have been developed to describe the response functions of photosynthesis to environmental factors.Here, we consider three main limiting factors imposed by radiation, temperature, and water on vegetation photosynthesis.The equations developed for the dual-source dualleaf (DSDL) model (Ding et al., 2014), Terrestrial Ecosystem Model (Raich et al., 1991), and Biome-BGC (BioGeochemical Cycles) models (Running et al., 2004) are used to account for the influences of radiation, temperature, and vapor pressure deficit (VPD), respectively: where k C and k Q are the stress coefficients of PAR absorbed by plant leaves; Q is the mean APAR for sunlit or shaded leaves per leaf hemi-surface area; T min , T opt , and T max are the minimum, optimum, and maximum temperature for photosynthetic activities, respectively; and VPD min and VPD max are the minimum and maximum vapor pressure deficit, respectively.In the DSDL model, k C and k Q are 500 and 150 W m −1 , respectively.T min , T opt , and T max are determined as 10, 28 and 48 • C for C4 crops (Kalfas et al., 2011), and here we slightly lower their values to 0, 25, and 45 • C, respectively, for C3 plants.VPD min and VPD max are 0.65 and 4.6 kPa for deciduous forests, respectively, in the Biome-BGC model (Heinsch et al., 2003).
Due to different PAR absorption by sunlit and shaded leaves, the stomatal conductance for sunlit and shaded leaves need to be calculated separately as where g cSun and g cShd are the stomatal conductance for sunlit and shaded leaves, respectively, and Q Sun and Q Shd are the mean PAR absorbed by sunlit and shaded leaves, respectively.
Given measured ambient CO 2 concentrations, the closure of the Eq. ( 47) now requires the quantity of intercellular CO 2 concentrations.Katul et al. (2000) compared eight models and concluded that all reproduced the measured carbon assimilation rates well.Here, we employ Leuning's method (Leuning, 1995) to estimate the ratio of intercellular to am- bient CO 2 concentrations as where VPD is the ambient vapor pressure deficit; VPD 0 is an empirical constant describing the species sensitivity to ambient vapor pressure deficit; is the leaf CO 2 compensation point; C a and C i are ambient and intercellular CO2 concentrations, respectively; and m L represents linear regression coefficients related to tree species.Calibrated values for model parameters are m L = 4.0, = 40 µmol mol −1 , and VPD 0 = 30 kPa, respectively (Katul et al., 2000).
Given modeled carbon assimilation rates at the leaf level, the total rate of carbon assimilation at the canopy level can be scaled up as where GPP is canopy gross primary production, A Sun and A Shd are leaf-level carbon assimilation rates for sunlit and shaded leaves, respectively, and LAI Sun and LAI Shd are the sunlit and shaded leaf area indexes.

Study materials and model parameterization
We studied two deciduous forest sites: Harvard Forest (US-Ha1) in Massachusetts and Bartlett Experimental Forest (US-Bar) in New Hampshire (Richardson et al., 2012).Basic information is briefly summarized in Table 1 for each site.Although plot layouts set up for the fieldwork did not match the exact footprints of flux towers (Yang et al., 2013), the measured tree structural attributes, such as tree density, are assumed to be representative of the two study sites.Flux towers measure energy and material fluxes between ecosystem and the atmosphere continuously (Baldocchi et al., 2001).Measured data are provided as standard Level 2 products in the AmeriFlux database (http://ameriflux.ornl.gov/).The time steps of available data are half-hourly for US-Bar and hourly for US-Ha1.The measurements we used include estimates of gross primary production (GPP) derived with the eddy covariance technique (Baldocchi, 2003), and meteorological variables such as shortwave solar radiation, temperature, vapor pressure deficit, and canopy-scale CO 2 concentration.Raw measurements of meteorological variables were used for analysis and missing values due to instrument malfunction or unsuitable micrometeorological conditions were screened.However, we obtained GPP estimates from AmeriFlux Level 4 products if they were not delivered in Level 2 products.Extraterrestrial solar radiation and solar zenith angle (i.e., the angle that the sun away from directly overhead) are calculated as a function of geolocation (i.e., latitude and longitude), the day of year (DOY), and solar time of the day (Allen et al., 1998).If diffuse radiation is missing from the measurements, we implement Muneer's method to partition global solar radiation into beam and diffuse components (Muneer, 2007): where f b is the proportion of beam radiation in global incoming radiation, and K t is the hourly clearness index.K t = I 0 /I e , where I 0 is global solar radiation on the canopy top and I e is the extraterrestrial solar radiation.We use typical parameter values from the literature for model parameterization.Because the spectral signatures of vegetation leaves and soil background differ in the spectral bands of PAR and near infrared (Table 2), we perform model simulations for these two discrete bands separately.Incident PAR is estimated to account for 47.5 % of incoming shortwave solar radiation, and the rest is attributed to the near-infrared band (Zhao et al., 2005).Maximum leaf stomatal conductance to H 2 O is estimated as 5.5 for US-Bar and 7.2 mm s −1 for US-Ha1 (Bonan, 2002;Ding et al., 2014), and they are translated to maximum leaf stomatal conductance to CO 2 assuming that the temperature is 20 • C and the atmospheric pressure is 101.32 kPa (Pearcy et al., 1989).Heights for canopy top (z 2 ) were measured to be 23.0 m for US-Ha1 and 19.0 m for US-Bar (Table 1), and heights for canopy bottom (z 1 ) were estimated as z 1 = 0.15 z z .Canopy structure in GORT is modeled with the ratios H /b = 2.0 and b/R = 3.0 (Strahler et al., 1999).Parameter values defined for canopy structure are somewhat arbitrary but are identical to our previous modeling efforts (Liu et al., 2008;Xin et al., 2012).The effects of tree structural parameters on model simulations are further explored in our study by varying their values.
Model validation for vegetation photosynthesis is performed with time series data for 8 successive days and for entire years.Based on AmeriFlux biological data, measured LAIs were 4.7 ± 0.   satellite-derived LAI from the MODIS (Moderate Resolution Imaging Spectroradiometer) products (Myneni et al., 2002).
The standard MODIS products (MOD15A2) provide 8-day LAI estimates at 1000 m spatial resolution, and we derived 8-day mean LAI for a 3 × 3 pixel window centered at each site.We screened cloudy observations based on the quality control data in MOD15A2 and applied double logistic equations to fit time series of cloud-free LAI observations (Li et al., 2014;Zhang et al., 2003).

Gap probability
The gap probabilities derived from the GORT model are shown in Fig. 2. As the solar zenith angle increases, more beams of sunlight are intercepted by leaves and tree crowns, resulting in decreased gap probabilities for both betweenand within-crown gaps.As LAI increases, within-crown gaps decrease but between-crown gaps remain the same.The physical explanation underlying is simple: tree leaves are clumped within each individual crown such that variations in LAI would not affect between-crown gaps, which are only a function of crown shape, canopy structure, and illumination geometry.Figure 3 further compares the gap probabilities modeled using GORT and Beer's law.For both models, gap probabilities decrease as solar zenith angle increases (Fig. 3a).Modeled gap probabilities are close when canopy LAI is low.However, at high LAI, the total gap derived from GORT is considerably greater than that modeled using Beer's law due to strong clumping effects.With an LAI of 4.0, the differences in gap probabilities are as much as 0.3 at the nadir and, in this case, more sunlight is allowed to be transmitted to the ground surface in GORT than in classic radiative transfer models.Modeled vertical structures of sunlight penetration are also shown to be different between GORT and Beer's law (Fig. 3b).The gap probability modeled using Beer's law decreases exponentially as canopy depth increases, whereas the decrease in the GORT-modeled gap probability follows an inverse sigmoidal curve.The reason behind this can be explained by the geometric factor: classic radiative transfer models assume that leaves are randomly distributed within the canopy layer, but the GORT model assumes that leaves are randomly distributed within individual crowns.Due to the ellipsoidal shape of tree crowns, there are simply more leaves in the canopy center than near the canopy top and canopy bottom, where the gap probability decreases more slowly.measured values (R 2 = 0.981), demonstrating the ability of GORT to model radiation absorption at the US-Bar site.

Model simulations over 8-day time periods
Time series of each component for modeling canopy photosynthesis are shown in Fig. 5. Given that total LAI remains the same over the course of several days, modeled sunlit and shaded LAI have little day-to-day variability and only vary as a function of solar zenith angle (Fig. 5a).As solar zenith angle decreases, sunlit LAI increases but shaded LAI decreases.Because sunlit leaves receive more illumination, they have less radiation limitations on photosynthesis than shaded leaves (Fig. 5b).Temperature limitation generally decreases from morning until noon, while VPD limitation increases.Although the chemical process of photosynthesis favors higher temperatures, leaf stomata tend to close to reduce water loss when atmospheric dryness is high (Bonan, 2002).Because short-term canopy CO 2 concentrations vary with winds and convection between the ecosystem and the atmosphere, the ambient CO 2 concentrations exhibit the greatest variation from day to day (Fig. 5b), so do the modeled differences between ambient and intercellular CO 2 concentrations.Figure 6 shows time series of measured and modeled GPP for two sites over 8 successive days.GPP estimates match flux tower measurements well in terms of the phase and amplitude.Daily peak GPP from tower measurements are over 30.0 µmol CO 2 m −2 s −1 for both sites.It is also evident that modeled results can capture some subtle variations in GPP at the hourly timescale.However, GPP estimates are slightly higher on DOY 242 but lower on DOY 243 for US-Ha1.Note that we used Muneer's method for estimating the diffuse radiation in US-Ha1 because measurements were not available.Considering uncertainties from the partition of global solar radiation, results for both sites perform well in general.
Figure 7 statistically compares measured and modeled GPP.Our model is able to explain 84.0 and 88.3 % of the GPP variances for the US-Bar and US-Ha1 sites, respectively.The regression lines are close to the 1 : 1 lines, and GPP is only slightly overestimated for US-Bar and underestimated for US-Ha1.The root mean squared errors (RM-SEs) are 3.71 and 3.08 µmol CO 2 m −2 s −1 for US-Bar and US-Ha1, respectively.The overall model performance is high considering that we did not attempt to perform model calibrations.

Model simulation over entire years
LAIs derived from satellite observations (Fig. 8) are used as inputs to model daily GPP over an entire year in addition to the 8-day model simulations.The double logistic fitting lines are shown to reduce noises in time series of MODIS LAI due to the effects of clouds and solar and viewing geometries.Fitted LAI time series are slightly higher from June to August and lower from September to December in 2006 at the US-Ha1 sites, but match with field measurements in general.The differences are likely to be introduced by mismatched observation footprints and uncertainties in satellite retrieval algorithms.The fitted time series of MODIS LAI are used for subsequent model simulations.
Figure 9 presents time series of measured and modeled GPP at the US-Bar site.Modeled results capture the trend and subtle variations of measured GPP on a daily basis.Most of the dips in the GPP time series occur on cloudy days when radiation is the main limiting factor for vegetation photosynthesis.GPP values at US-Bar are slightly overestimated from DOY 100 to 150 in 2004 possibly due to overestimation of the LAI.Statistically, modeled results can explain 79.5, 89.7, and 89.3 % of the variance in daily GPP for the years 2004, 2005, and 2006, respectively (Fig. 10).Regression slopes are close to the 1 : 1 lines except in the year 2004 due to overestimated GPP in the early growing season.The RMSEs are 1.64, 1.31, and 1.56 gC m −2 day −1 for 2004, 2005, and 2006, respectively.
Because measurements of atmospheric CO 2 concentrations within the canopy are largely unavailable for US-Ha1 (only approximately 41.4 % of the measurements are valid for use), we do not aggregate hourly results to daily sums but perform regression analysis using all available hourly data in Fig. 11.For the US-Bar site, the R 2 value is 0.801 and the RMSE value is 4.31 µmol CO 2 m −2 s −1 .For the US-Ha1 site, the correlation between modeled and measured GPP is strong with an R 2 value of 0.777 and an RMSE value of 6.49 µmol CO 2 m −2 s −1 .There were slight GPP underestimates when measured GPP values were high at the US-Ha1 site, possibly due to empirical functions that we used in modeling diffuse radiation and leaf photosynthesis.Table 3 lists major statistical results for our model performance, as evaluated using all available hourly data at both sites.The model performance is consistent through time and is comparable to the simulation of 8-day data (Fig. 7), despite the fact that satellite-derived LAI instead of field measurements were used for yearly simulation.

Influence of CO 2 concentration on canopy photosynthesis
One important question is whether it is necessary to link radiative transfer with leaf stomatal conductance for modeling photosynthesis, since some biogeochemical models such as production efficiency models simply assume that vegewww.biogeosciences.net/12/3447/2015/Biogeosciences, 12, 3447-3467, 2015 tation GPP/NPP is linearly related to canopy radiation absorption (Xin et al., 2013).To understand the performance of production efficiency models, we conduct linear regressions between modeled APAR and measured GPP as shown in Fig. 12.Indeed, canopy APAR is positively related to flux tower GPP and explains 70.3 % of its variance.The R value increases slightly to 0.710 after accounting for the influences of temperature and vapor pressure.The model performance here is comparable to results from other studies that evaluate production efficiency models (Chen et al., 2011;Sjöström et al., 2013;Xin et al., 2015).However, there are strong partial correlations between canopy CO 2 concentrations and GPP even after accounting for radiation absorption.Figure 13a shows the residual plot of GPP versus ambient CO 2 concentrations when controlling on APAR.The slope is negative because the am- The data clearly allow for rejection of the null hypothesis that ambient CO 2 concentration has no effects on canopy photosynthesis.This relationship holds even after considering the factors of temperature and vapor pressure deficit (Fig. 13b).
We therefore conclude that accounting for the influence of ambient CO 2 concentrations is essential for modeling daytime GPP at the half-hourly timescale.

Clumping effects in the GORT model
The clumping effects of leaves modeled using GORT influence canopy radiative transfer processes and are worthy of further examination.Chen et al. (1997) demonstrated that the net effects of leaf clumping could be modeled by introducing a clumping index.We derive the clumping index by inverting their functions (Zhao et al., 2011) as follows: where is the clumping index, P gap is the gap probability modeled using GORT, P Beer = exp(−k b LAI) is the gap probability modeled using Beer's law, k b is the extinction coefficient, and LAI is the leaf area index.The behavior of the derived clumping index shown in Fig. 14 is intuitively interpretable.Leaves are more clumped when LAI is larger given constant tree structures.However, when LAI is constant but tree density increases, leaves are distributed in a larger three-dimensional space, resulting in an increased clumping index.Similarly, if the H /b ratio or b/R ratio decreases while other parameters are unchanged, the total crown volume increases and leaves are less clumped.The sensitivity of the clumping index to the illumination zenith angle varies when using different parameter sets.Our simulated results are in line with the measured and modeled results in previous studies (Leblanc and Chen, 2001;Leblanc et al., 2002): the clumping indexes are insensitive to zenith angles in some forest stands and increase with zenith angles in others.We do not attempt to derive clumping indexes at solar zenith angles greater than 85 • when gap fractions typically approach zeros.These results imply that tree structure strongly influences radiation absorption and photosynthesis of canopies.

Assumptions and future improvements
It is also necessary to review our model assumptions and identify possible avenues for future improvements.First, we assume a spherical leaf angle distribution in the model simulations.However, most deciduous forests have semihorizontal leaf orientation (Bonan, 2002) and an assumption of planophile or plagiophile LAD is likely to be more appropriate for temperate and boreal broadleaf forests (Pisek et al., 2013).Because LAD influences the proportions of sunlit and shaded leaf areas, the way in which modeled canopy GPP varies with LAD requires further exploration.Second, the substrate under the canopy layer is assumed to be a Lambertian surface.Field studies have observed the effects of bidirectional reflectance distribution function (BRDF) for soils (Liang and Townshend, 1996;Wang et al., 2010), and coupled soil and vegetation models (Ni and Li, 2000;Verhoef and Bach, 2007) should be tested to understand the effects of soil BRDF on canopy photosynthesis.Third, we assume maximum constant leaf stomatal conductance over the growing season.It is worth examining how optimal leaf stomatal conductance may evolve with leaf development stages and long-term environmental changes (Keenan et al., 2013;Lammertsma et al., 2011).Fourth, we use ellipsoids to describe tree crown shapes for deciduous broadleaf forests.Because many evergreen needleleaf forests have conical crowns, future applications to areas with conifer forests may require different treatment for crown shapes in the models.Fifth, multi-story vegetation canopy such as overstory and understory are common in forest ecosystems, and it may be neces-  sary to improve the current model by considering multi-layer vegetation canopies in future studies.Finally, our linkage between radiative transfer and biochemical processes is still empirical.We may need to test other mechanisms, for example the biochemical model based on the enzyme kinetics of rubisco (ribulose bisphosphate carboxylase/oxygenase) and the regeneration of RuBP (ribulose-1, 5-bisphosphate) in response to light absorption (Farquhar and Sharkey, 1982), in future studies.

Conclusions
We propose and validate a new model that links GORT with biochemical processes for modeling canopy photosynthesis.Several main conclusions can be drawn from this study.First, the radiative transfer process within the canopy is one of the key factors in modeling vegetation photosyn-thesis, and our proposed model simulates canopy photosynthesis well.Modeled GPP robustly explained approximately 80 % or more variance in GPP measurements at both halfhourly and daily timescales.Second, tree structures influence canopy gap probabilities and vegetation photosynthesis.Leaf clumping could vary as a function of tree density, canopy depth, and crown shapes and affect canopy sunlight interception.Finally, ambient CO 2 concentrations influence vegetation photosynthesis activities and should be included in biogeochemical models.

Figure 1 .
Figure 1.A scheme of the canopy structure in the Geometric Optical Radiative Transfer model as modified from Ni (1998).
2 on DOY 211 in 2004 at the US-Bar site and 4.84 ± 0.78 on DOY 234 in 2006 at the US-Ha1 site.Because field-measured LAI data were insufficient to support model simulation for an entire calendar year, we obtained www

Figure 2 .
Figure 2. Canopy gap probabilities modeled using GORT with varied leaf area index.The total gaps are between-crown gaps plus within-crown gaps.Tree structure parameters for the US-Bar site are used in model simulation.

Figure 3 .
Figure 3. Comparisons between canopy gap probabilities modeled using GORT and Beer's law as a function of (a) solar zenith angle and (b) canopy depth.The canopy depth is defined as the distance from canopy top to a canopy height (h).Tree structure parameters for the US-Bar site are used in GORT simulation.

Figure 4 .
Figure 4. Measured and modeled components of radiation in 8 successive days are shown for (a) the partition of global solar radiation, (b) surface radiation balance, (c) modeled and measured diffuse radiation, and (d) modeled and measured net radiation.Extraterrestrial radiation is derived following methods outlined in Allen et al. (1998).Muneer's method is applied to model diffuse radiation.The GORT model is applied to model net radiation.Data are shown from DOY 217 to 224 in 2004 for the US-Bar site.

Figure 4 Figure 5 .
Figure4shows each component of the radiation regime at the US-Bar site.The diffuse radiation modeled using Muneer's

Figure 6 .
Figure 6.Time series of modeled and measured GPP for 8 consecutive days at the sites (a) US-Bar and (b) US-Ha1.Data are halfhourly at the US-Bar site and hourly at the US-Ha1 site.Data are shown from DOY 217 to 224 in 2004 for US-Bar, and from DOY 241 to 224 in 2006 for US-Ha1.Negative GPP measurements are set to zero.Missing points in modeled GPP at the US-Ha1 site are due to missing measurements of canopy CO 2 concentrations or other meteorological variables.

Figure 7 .
Figure 7. Regressions between modeled and measured GPP for 8 consecutive days at the sites (a) US-Bar and (b) US-Ha1.Data are from DOY 217 to 224 in 2004 for US-Bar and from DOY 241 to 224 in 2006 for US-Ha1.Only data during the photosynthetically active period (flux tower GPP > 0.5 µmol CO 2 m −2 s −1 ) are included in the regression.The solid lines denote the 1 : 1 lines, and the dashed lines denote the regression lines.

Figure 8 .
Figure 8. Comparisons of field-measured and satellite-derived LAIs for the sites (a) US-Bar in 2004 and (b) US-Ha1 in 2006.The solid grey lines denote MODIS LAI as obtained from standard MODIS FPAR/LAI products (MOD15A2).The solid black lines denote double logistic fitting lines that are applied to MODIS LAI.The solid points denote the measured LAI as obtained from biological data sets from the AmeriFlux website.

Figure 9 .
Figure 9.Time series of modeled and measured daily GPP shown for (a) 2004, (b) 2005, and (c) 2006 at the US-Bar site.Model simulation is performed at a half-hourly time step.Measured and modeled half-hourly GPP are aggregated to generate daily time series with units converted from µmol CO 2 m −2 s −1 to gC m −2 day −1 .Occasional negative GPP measurements are set to zeros.Missing points in modeled GPP time series are due to missing measurements of meteorological variables during the daytime photosynthetically active period (flux tower GPP > 0.5 µmol CO 2 m −2 s −1 ).

Figure 10 .
Figure 10.Regressions between modeled and measured daily GPP shown for (a) 2004, (b) 2005, and (c) 2006 at the US-Bar site.Only data during the photosynthetically active period (flux tower GPP > 0.5 g C m −2 day −1 ) are included in the regressions.The solid lines denote the 1 : 1 lines, and the dashed lines denote the regression lines.

Figure 11 .
Figure 11.Regressions between modeled and measured GPP for all available hourly data at the sites of (a) US-Bar and (b) US-Ha1 in 2006.Only data from the photosynthetically active period are included in the regression.The solid lines denote the 1 : 1 line, and the dashed lines denote the regression line.

Figure 12 .
Figure 12.Regressions between modeled APAR and measured GPP.Half-hourly data are shown from DOY 217 to 224 in 2004 for US-Bar.The influences of temperature and vapor pressure deficit are modeled based on Equations (51) and (52).Only data during the photosynthetically active period are included in the regression.The dashed lines denote the regression lines.

Figure 13 .
Figure 13.Residual plots are shown for (a) the partial correction between GPP and ambient CO 2 concentration (C a ) while controlling for the variable of APAR and (b) the partial correction between GPP and C a −C i while controlling for the variable of APAR ×f (T ) × f (VPD).

Figure 14 .
Figure 14.Derived clumping index as a function of solar zenith angle for varied canopy parameters.Tree parameters for US-Bar are used for GORT simulations.The default simulation is for a canopy composed of H /b = 2.0, b/R = 3.0, λ = 1432 trees ha −1 , and LAI = 2.0, and labeled curves are for the same case with only the labeled parameters varied.

Table 1 .
Site information as obtained from the AmeriFlux website unless notified.

Table 2 .
The spectral signature of leaf and soil background.

Table 3 .
The model performance at two study sites as evaluated using hourly data.Units for RMSE and mean bias error (Bias) are in micromoles of CO 2 per square meter per second (µmol CO 2 m −2 s −1 ).CO 2 concentration, as regulated by vegetation photosynthesis and respiration activities, is normally high during the nighttime but low during the daytime.The correlation coefficient is only −0.279, but it is statistically significant (p value < 0.001) under a one-tailed partial correlation test. bient

Table A1 .
transmittance for semi-infinite homogeneous canopies ρ ff (h) hemispherical-hemispherical reflectance for homogeneous canopies with finite thickness ρ ff (h) directional-hemispherical reflectance for homogeneous canopies with finite thickness t ff (h) hemispherical-hemispherical transmittance for homogeneous canopies with finite thickness t df (h, θ i ) directional-hemispherical transmittance for homogeneous canopies with finite thickness ρ ff (h) hemispherical-hemispherical reflectance for discontinuous canopies ρ df (h, θ i ) directional-hemispherical reflectance for discontinuous canopies t ff (h) hemispherical-hemispherical transmittance for discontinuous canopies t df (h, θ i ) directional-hemispherical transmittance for discontinuous canopies δLAI (h) leaf area index within a thin layer δh at height h LAI total leaf area index of the canopy LAI Sun (θ i ) sunlit leaf area index given a solar illumination angle θ i LAI Shd (θ i ) shaded leaf area index given a solar illumination angle θ i LAI * Sun (θ i ) sunlit leaf area for homogeneous canopies given a solar illumination angle θ i Continued.Symbols Definition f b the fraction of incident beam radiation in total or global incoming solar radiation I b (h, θ i ) unintercepted beam fluxes at canopy height h given a solar illumination angle θ i I d (h)unintercepted diffuse fluxes at canopy height h I bt (h, θ i ) unintercepted and down-scattered beam fluxes I dt (h) cmax maximum leaf stomatal conductance when environmental factors do not limit carbon uptake f (x i ) scalars that account for the influences of environmental stresses on leaf stomatal conductance f (Q) scalars that account for the influences of solar radiation on leaf stomatal conductance f (T ) scalars that account for the influences of temperature on leaf stomatal conductance f (VPD) scalars that account for the influences of vapor pressure deficit on leaf stomatal conductance k C stress coefficients of PAR absorbed by plant leaves for the temperature scalar k Q stress coefficients of PAR absorbed by plant leaves for the temperature scalar Beer gap probability for beam light passing through the canopy as modeled using Beer's law www.biogeosciences.net/12/3447/2015/Biogeosciences, 12, 3447-3467, 2015

Table A2 .
Values for model parameters.