Using Stochastic Ray Tracing to Simulate a Dense Time Series of Gross Primary Productivity

Eddy-covariance carbon dioxide flux measurement is an established method to estimate primary productivity at the forest stand level (typically 10 ha). To validate eddy-covariance estimates, researchers rely on extensive time-series analysis and an assessment of flux contributions made by various ecosystem components at spatial scales much finer than the eddy-covariance footprint. Scaling these contributions to the stand level requires a consideration of the heterogeneity in the canopy radiation field. This paper presents a stochastic ray tracing approach to predict the probabilities of light absorption from over a thousand hemispherical directions by thousands of individual scene elements. Once a look-up table of absorption probabilities is computed, dynamic illumination conditions can be simulated in a computationally realistic time, from which stand-level gross primary productivity can be obtained by integrating photosynthetic assimilation over the scene. We demonstrate the method by inverting a leaf-level photosynthesis model with eddy-covariance and meteorological data. Optimized leaf photosynthesis parameters and canopy structure were able to explain 75% of variation in eddy-covariance gross primary productivity estimates, and commonly used parameters, including photosynthetic capacity and quantum yield, fell within reported ranges. Remaining challenges are discussed including the need to address the distribution of radiation within shoots and needles.


Introduction
Monitoring forest productivity not only widens our ecological understanding, but is increasingly used as a diagnostic tool for evaluating silvicultural practices and testing compliance with political regulations regarding sustainable resource use.From an ecological point of view, information about forest productivity is often indicative of biodiversity richness [1] and can be used to explain how anthropogenic influences and climatic changes affect forest health and functioning [2].A key variable is the gross primary productivity (P g ) that defines the uptake of carbon from the atmosphere before any losses due to respiration are subtracted.Estimates of P g are based on the amount of photosynthetically active radiation that is absorbed by the canopy (E ac ) and its pathways in photosynthesis [3][4][5].The modelling of P g is challenging, however, due to the size and variability of the global forest resources, the multitude of environmental drivers, and the lack of control over these drivers in the field.As a result, a wide range of models exists that operate at a wide range of temporal and spatial scales.To validate these models, a hierarchical strategy is generally adopted where P g models at broad scales are calibrated and validated using models and sampling techniques at finer scales.
At broad spatial scales, such as regional to global, forest productivity models represent the canopy either as a single big leaf that extends across the landscape [3,4,6] or as a stack of foliage layers that provide for the modelling of vertical transfer of visible and thermal radiation [7,8].At these model scales, the net lateral transport of energy is assumed to be zero and the models are often parameterized using remote sensing estimates of incident short-and long-wave radiation and stand-level estimates of canopy structure, site nutrition, and hydrology [6,[9][10][11].Paramount to the calibration and validation of these models is the eddy-covariance (EC) method that provides estimates of P g at stand level (P gs ) from high-frequency measurements of the vertical wind velocity, and carbon dioxide (CO 2 ) and water vapour (H 2 O) mixing ratios [12].P gs is obtained by subtracting daytime EC-measured net ecosystem exchange of CO 2 (NEE) from daytime ecosystem respiration (R e ) usually estimated using the relationship between nighttime EC-measured NEE and soil or air temperature.Currently, a global network of over 500 EC installations provides biome-specific, long-term information about variability in P gs within footprints of around 10 ha.
For homogeneous forest stands on horizontal terrain, the error bound on NEE obtained with the EC method is less than 50 g¨C¨m 2 ¨year ´1 with long-term measurements and independent flux measurements of individual ecosystem components adding additional confidence [13].For example, measurement of CO 2 exchange can be undertaken at scales ranging from leaf to tree level (using chambers or cuvettes) and can be used in support of EC-data acquisition.At these fine scales, the heterogeneity of the canopy structure requires a consideration of vertical as well as lateral transport of energy [14].To compute the three-dimensional (3D) canopy radiation field, an analytic approach is widely used that represents the stand as a series of conical and ellipsoidal crown envelopes for which smooth, radial transitions in leaf density can be specified (e.g., [5,[15][16][17].More recent approaches have adopted stochastic ray tracing to simulate the scattering of light within increasingly heterogeneous plant canopies, based on first principles [18][19][20][21].In ray tracing, the scattering and absorption of individual photons is simulated by computing intersections between light rays and the scene and by sampling light propagation from material-specific bi-directional reflectance (transmittance) distribution functions (BRDF and BTDF, respectively).Ray tracing enables the simulation of absorbed radiation for individual canopy elements (E al ) and holds significant potential to scale P g estimates from the stand level to finer scales where experiments can be increasingly controlled.Virtual 3D canopies can be generated using tree-growing algorithms that are based on a recursive application of a set of growing rules [22][23][24].To limit data dimensionality, a clump of leaves may be represented by its convex hull [24] or as a planar polygon or plate-like shoot model that exhibits scattering properties similar to the assembly of individual leaves or needles [25,26].Recent studies have demonstrated the combined use of tree-growing algorithms and time-of-flight laser scanning technology [27,28] for the reconstruction of virtual forest canopies [29][30][31][32][33].
Ray tracing simulations are computationally expensive and require the use of accelerating structures that group scene elements into hierarchical bounding volumes or regular grids to speed up ray-intersection tests [34,35].These accelerating structures provide performance gains for fixed illumination conditions; however, they do not provide additional gains for cases where illumination conditions are dynamic, i.e., when the distribution of radiation across leaves is to be studied over time and for different sun angles and cloudiness conditions.As a result, the application of ray tracing for scaling between P gl and P gs estimates over extensive time series remains challenging.
The work presented herein relies on the modelling of time series of canopy radiation transfer using stochastic ray tracing.To speed-up ray tracing simulations under dynamic lighting conditions, we compute a look-up table designed to hold the probabilities of absorbing radiation from over a thousand hemispherical directions for thousands of scene elements.Once this look-up table is populated, the absorption of photosynthetically active radiation is computed for each element, with only array multiplications, thereby significantly reducing computing costs.To improve computation time, the virtual forest scene is modelled at a reduced resolution of 50,000 elements, representing soil, bark and foliage.We demonstrate the inversion of this model and estimate empirical leaf-level photosynthesis parameters from EC-estimates of P gs and meteorological measurements.We conclude by discussing the limitations of this reduced scene fidelity and challenges and opportunities ahead for the proposed technique.

Study Area
The study area is located on Vancouver Island, BC, Canada, near the city of Campbell River, in the Oyster River area about 16 km from the east coast of the island (49 ˝52'7.8"N,125 ˝20'6.3"W,350 m above mean sea level).The area is a plantation forest consisting of 80% Douglas-fir (Pseudotsuga menziesii spp.menziesii (Mirb.)Franco), 17% western redcedar (Thuja plicata Donn.ex D. Don), and 3% western hemlock (Tsuga heterophylla (Raf.)Sarg) [36][37][38] with main understorey species, salal (Gaultheria shallon Pursh.), dull Oregon grape (Mahonia nervosa (Pursh.)Nutt.), vanilla-leaf deer foot (Achlys triphylla (Smith) DC), and various ferns and mosses.The site is located on a northeast facing 5 ˝-10 ˝slope and has a stand density of 1100 stems ha ´1 with tree heights ranging between 30 and 35 m [39], a mean diameter at 1.3 m height of 31 cm [40] and a mean annual increment of 12.1 m 3 ¨ha ´1 [38].The predominant age of the trees was 60 years old in 2009 and the leaf area index was 7.3 m 2 ¨m´2 [41].The soil is classified as a humo-ferric podzol and has a gravelly sandy loam texture and a surface organic layer ranging in depth between 1 and 10 cm [38].The region belongs to the dry maritime Coastal Western Hemlock biogeoclimatic subzone (CWHxm) and has a mean annual precipitation and temperature of 1470 mm and 8.6 ˝C, respectively [40].Approximately 110 ha of the study area, extending in the prevailing wind directions (NE and SW), was aerially fertilized with urea at 200 kg¨N¨ha ´1 on 13 January 2007 and simulation of the footprint climatology has demonstrated that this extent contributes to over 80% of the eddy-covariance flux measured at the flux tower [40,42].

Data
One 30 ˆ30 m plot was selected based on its representativeness of the stand conditions [43] and its close proximity to the eddy-covariance fluxtower.In 2008, ground-based laser scanning data were acquired with the EVI (CSIRO) full-waveform laser scanner that emits roughly half a million laser pulses in an upward hemispherical direction and records the reflectance of this laser light at a sampling rate of 2 Gs/s so that the distance to the scattering surfaces can be determined at centimeter accuracy, see e.g., Strahler et al. [44] for additional details.By emitting pulses over a more than hemispherical field of view (i.e., 130 ˆ360), millions of returns are recorded per scan and a total of five scans were taken: one in the centre and one per plot corner, to provide for a comprehensive coverage of plot-level forest structural information relating to stem shape and leaf area distribution.Half-hourly NEE along with meteorological records including incident diffuse and total photosynthetically active radiation (E ic ), relative humidity and temperature were selected for the year 2009.NEE was calculated as the sum of above-canopy EC-measured CO 2 flux and the estimated rate of change in CO 2 storage in the air column beneath the EC sensors [38].Diffuse and total E ic were measured using a BF2 sunshine sensor (Delta-T Devices, Cambridge, UK) mounted on top of the EC-flux tower (at 45 m).Air temperature and relative humidity were measured at the flux tower at the 27-m height.P gs was calculated as daytime R e minus daytime NEE.The former was estimated using the exponential relationship between measured nighttime NEE and soil temperature at the 2-cm depth [37].

Scene Reconstruction
A 3D virtual forest scene was reconstructed from the ground-based laser scanning data using a modeling approach described by [33].Briefly, single scans are projected onto a panoramic plane where individual pixels show range (m) to target.Individual tree stems are detected from the panoramic projection using image segmentation and filtering criteria.Image segmentation is used to detect surface edges and continuous surfaces that show sharp and gradual changes in range, respectively.Continuous surfaces relate to ground, foliage and stem objects and hence filtering criteria, including auto-correlation along surface edges and segment lengths are used to separate stem segments from other targets in the scans.Three-dimensional stem objects were created using the range information and diameter estimates derived from the laser scanning data.All scans and stem objects were merged to a common coordinate system that was clipped to the plot boundaries (i.e., 30 m ˆ30 m).Overlapping stem objects were merged into a single object that was assigned the average stem attributes obtained from the individual scans.From the co-registered point cloud, a canopy height model (CHM), a digital elevation model (DEM), and a delineation of individual crowns [45] were obtained.Gaps due to data occlusion along the stem segments were filled using linear interpolation and the stem tips were extrapolated by tracing the local point cloud geometry [46].Using the stem information, individual tree growing spaces were derived based on 2D Voronoi tesselation for horizontal slices across a regular vertical interval.A template tree crown was constructed using Arbaro [23] (Figure 1): computer animation software for the realistic modelling of crown architecture based on simple geometric parameters such as branch lengths and insertion angles of child-branches in relation to their parent branches.Scaled copies of the crown template were then fitted to the virtual scenes.In contrast to skeletal tree structures produced by Arbaro that contain a high level of structural detail (Figure 2a), the template tree used in this study had a much coarser resolution to reduce computation time.Recent studies have explored efficient means to abstract shoots and whorls using a series of projection planes (e.g., [25,26]).The template tree used in this study was constructed using a layered series of planar polygons producing a much reduced scene complexity (Figure 2b).This abstraction was chosen to demonstrate the proposed concept, acknowledging that a comprehensive treatment of needle-level details at plot scale requires dedicated computing facilities and more advanced acceleration techniques.The chosen level of abstraction foremost affects transmissions from oblique incidence angles, as well as within-shoot distributions of E al and re-collision and escape probabilities (e.g., [26]).Effects on the latter are considered small, however, in the PAR band.For each first-order branch, a planar polygon was fit to its principle stem and all of its sub-stems.This fit was computed as a two-dimensional triangulation of vertices (i.e., height-coordinates were not used in triangulating).The structural fidelity was chosen to resemble the clumping of foliage into crowns, but did not include clumping of needles into shoots or clumping of shoots around whorls [47,48].To provide for the simulation of un-collided transmission within whorls, gap fractions (GF K ) were assigned to the facets of the planar polygons.Whorl-level gap fraction was estimated as the fraction of visible sky (15%) in digital photographs taken underneath isolated Douglas-fir branches [33].To assess model sensitivity to GF K , three scenarios were considered throughout the remaining analysis, for gap fractions of 10%, 15%, and 20%, respectively.Un-collided transmission τ u amounts to the gap fraction when the incidence angle is perpendicular and drops with the cosine of the incidence angle.Finally, the combined area of all planar polygons was chosen to match 1{ p1 ´0.15q of stand LAI.The 3D scene was stored in the OBJ-file format (formerly Wavefront technologies, now Alias Systems Corp., Toronto, ON, Canada) as a sequence of individual, connected triangles, comprising a total of 102 trees.The integrity of the virtual scene was validated using full-waveform derived canopy gap probabilities [49] and gap probabilities derived through simulation (see [33], for further details).Validation of the ray tracer was achieved by means of conservation-of-energy tests, introspection of traced paths using 3D interactive visualization, and four RAMI baseline test cases relating to foliar absorptance and canopy transmission [33].

Canopy Radiation Modeling
The virtual scene can be used in a ray tracer that simulates the propagation of individual photons and their interactions with the scene.Photons travel along lines that may intersect with surface elements.At the location of the nearest intersection, a fate is computed as to whether the photon is absorbed, scattered, or penetrates the element freely without collision.A constant Lambertian scattering probability was used across all scene elements representing foliage, , ( = 0.05, , = 0.02) and a constant Lambertian reflectance for all soil and bark elements, and respectively ( = = 0.1).For a Lambertian surface, the direction of scattering (i.e., reflection or transmission) is independent from the photon incidence angle and follows a cosine-weighted hemispherical probability distribution where the reflected intensity drops with the cosine angle [35,50].Earlier work [18,51] showed the use of ray tracers for the simulation of transmittance, reflectance, and absorptance by virtual plant canopies.In forward ray tracing simulations, photons are emitted from a horizontal reference plane that is located just above the 3D scene and interactions with the scene are computed for fixed illumination and viewing conditions.The requirement to simulate a large number of rays to

Canopy Radiation Modeling
The virtual scene can be used in a ray tracer that simulates the propagation of individual photons and their interactions with the scene.Photons travel along lines that may intersect with surface elements.At the location of the nearest intersection, a fate is computed as to whether the photon is absorbed, scattered, or penetrates the element freely without collision.A constant Lambertian scattering probability was used across all scene elements representing foliage, , ( = 0.05, , = 0.02) and a constant Lambertian reflectance for all soil and bark elements, and respectively ( = = 0.1).For a Lambertian surface, the direction of scattering (i.e., reflection or transmission) is independent from the photon incidence angle and follows a cosine-weighted hemispherical probability distribution where the reflected intensity drops with the cosine angle [35,50].Earlier work [18,51] showed the use of ray tracers for the simulation of transmittance, reflectance, and absorptance by virtual plant canopies.In forward ray tracing simulations, photons are emitted from a horizontal reference plane that is located just above the 3D scene and interactions with the scene are computed for fixed illumination and viewing conditions.The requirement to simulate a large number of rays to

Canopy Radiation Modeling
The virtual scene can be used in a ray tracer that simulates the propagation of individual photons and their interactions with the scene.Photons travel along lines that may intersect with surface elements.At the location of the nearest intersection, a fate is computed as to whether the photon is absorbed, scattered, or penetrates the element freely without collision.A constant Lambertian scattering probability was used across all scene elements representing foliage, p l , `ρl " 0.05, τ c,l " 0.02 ˘and a constant Lambertian reflectance for all soil and bark elements, p s and p b respectively pρ s " ρ b " 0.1q.For a Lambertian surface, the direction of scattering (i.e., reflection or transmission) is independent from the photon incidence angle and follows a cosine-weighted hemispherical probability distribution where the reflected intensity drops with the cosine angle [35,50].Earlier work [18,51] showed the use of ray tracers for the simulation of transmittance, reflectance, and absorptance by virtual plant canopies.In forward ray tracing simulations, photons are emitted from a horizontal reference plane that is located just above the 3D scene and interactions with the scene are computed for fixed illumination and viewing conditions.The requirement to simulate a large number of rays to obtain a realistic distribution of canopy radiation makes the use of ray tracers, challenging when illumination conditions are dynamic and updates to the radiation budgets are needed in real-time [52].The current method populates a look-up table (LUT) of which the rows correspond to scene elements ppq and columns to sky directions pφ i , θ i q.Each table element corresponds to the probability (P abs ) of absorbing photon flux (µmol¨s ´1) that passes the reference plane from direction pφ i , θ i q by canopy element p, so that the sum of probabilities across rows corresponds to canopy absorptance from direction pφ i , θ i q.Using this LUT, the irradiance that is absorbed by a scene element can be estimated by integrating radiance over directions pφ i , θ i q weighted by P abs :

E al ppq "
A re f A p ¨ P abs pp, φ i , θ i q ¨Lic pφ i , θ i q ¨cos pθ i q ¨sin pθ i q ¨dθ i ¨dφ i In Equation ( 1), A re f ¨ L ic pφ i , θ i q ¨cos pθ i q ¨sin pθ i q ¨dθ i ¨dφ i expresses the flux of photons received over a projected solid angle that passes through the reference plane that is located just above the highest scene element.Inclusion of P abs in the integral determines the fraction of photons passing the reference plane from the projected solid angle Ω `" cosθ i ¨sinθ i ¨dθ i ¨dφ i ˘that are absorbed by scene element p and accounts for the orientation of the elements as well as shading, and multi-scattering between the elements.P abs is computed as the fraction of photons from direction pφ i , θ i q that is absorbed by element p. Finally, division by A p normalizes this absorbed quantity to µmol photons m ´2¨s ´1.
To improve computing time, the sky was sampled using a Capitulum sampling distribution that owes its name to the resemblance to the distribution of seeds in a sunflower head.Samples were distributed about the zenith so that the distance to that axis decreased exponentially with the sample's rank and the angular distance between samples was approximately a constant of 2π{n (e.g., [53,54]).The hemisphere was covered with 1366 samples, so that the average solid angle was 0.0046 mrad.To avoid edge effects for the small study plot, cyclic boundaries were applied so that a ray from direction pφ i , θ i q escaping the scene at location px, y, zq was made incident from the opposite side of the scene.
Half-hourly radiances L (µmol¨photons¨m ´2¨s ´1¨sr ´1) were computed from diffuse and total E ic .Diffuse radiance was assumed to be homogeneous across the hemisphere, and direct radiance and sun position were computed using Ephem, a library for astronomical computations available to the Python programming language.Direct radiance was then assigned to the hemispherical sample that was located closest to the sun's position, determined using a KD-tree.Absorbed radiation at the canopy element level (µmol¨photons¨m ´2¨s ´1) was computed for every hour between 5 a.m. and 9 p.m. PST, for the majority of the growing season between 1 May and 17 September 2009.

Photosynthesis Modeling
A leaf-level photosynthesis model was developed based on a non-rectangular hyperbolic relationship between P gl and E al , following [55].This model is parameterized by photosynthetic capacity (P max ) that defines the asymptote, the quantum yield (α) that defines the initial slope when photosynthesis is electron transport limited, and the duration of the initial linear response (χ) that is larger for cells than for leaves and larger for leaves than for canopies (e.g., [55]): Down-regulation of photosynthesis is expected when the leaf is light saturated and when photosynthesis is limited, for example, by water availability or when the carboxylation rate of the Calvin cycle due to temperature or nutrient supplies is suboptimal.This effect of down-regulation was added to the photosynthesis model using air temperature (T) and relative humidity (h) modifiers that, for this purpose, were modelled as symmetric bell-shaped functions, acknowledging that investigations into asymmetry and alternative function shapes are outside the scope of this paper: Remote Sens. 2015, 7, 17272-17290 where T opt is the optimal temperature for photosynthesis, and T lag and h lag are the temporally lagged temperature and relative humidity conditions that were derived following [56].For temperature, this equation has the form: T lag ptq " p1 ´γT q ¨T ptq `γT ¨Tlag pt ´1q where T ptq represents the temperature at time t, and γ T defines the lag (γ T PR , 0ď γ T ď 1).The relative humidity modifier is analogous to that for temperature.Two additional modifier functions were tested to include a seasonal adjustment of photosynthetic capacity that corresponds to different phenological stages of tree growth and the effects of acclimation of the photosynthetic capacity to prevalent PAR [21,57]: k rub defines the decrease of nitrogen content with canopy depth and LAI is the cumulative leaf area index at facet p and is zero at the top of the canopy, D is growing degree days and was for this study computed as the mean daytime temperature minus 10 ˝C accumulated over days, starting from 1 May 2009.The down-regulation of P max and α was modeled as the product of an optimal value for conditions when photosynthesis is not limited and modifier functions: The resulting leaf model has at most eleven parameters (χ, P opt , α opt , T opt , D opt , β h , β T , β D , γ h , γ T , k rub ).Stand-level P gs for a specific moment in time t was computed through integration over all canopy elements n:

Model Inversion
In forward mode, stand-level P gs is obtained by integrating leaf-level assimilation over the canopy.In inverse mode, the leaf model parameters are inferred using the EC-estimates of P gs , T, h, and simulated E al .Model inversion is achieved numerically by minimizing the cost function: argmin Θ ÿ t " P gs, sim pΘ; E ic , T, h; tq ´Pgs, EC ptq ‰ 2 (10) where P gs,sim and P gs,EC are the simulated and EC-derived estimates of P gs obtained from Equation (9) and methods described in Section 2.2, respectively.A number of optimization techniques exist [58] and in this study, the Levenberg-Marquardt method was applied that minimizes the cost function iteratively and requires initial estimates of the parameter set Θ used to define leaf-level photosynthetic assimilation.To study the effects of different down-regulation functions, four cases were considered (Table 1): The first case ignores down-regulation and uses fixed values for P max and α across all photosynthetic scene elements, p l .The second case investigates down-regulation of photosynthetic capacity and quantum yield based on temperature and relative humidity conditions.The third case includes the seasonal adjustment of photosynthetic capacity and the fourth case includes the acclimation of foliage elements to sun and shade conditions, following [57].To establish convergence, χ was constrained to 0.9, which is slightly higher than 0.8 used by [17] yet within the range specified by [59] for leaf-level photosynthesis (i.e., 0.5 to 0.9).In addition, some parameter values used in cases three and four were fixed to estimates obtained from cases two and three, respectively.Temperature parameters (i.e., β T , γ T , T opt ) were fixed in case three and both temperature and relative humidity parameters (i.e., β T , γ T , T opt , β h , γ h ) in case four.Initial values for free parameters were 10 µmol¨C¨m ´2¨s ´1, 0.1 µmol¨C¨µmol ´1 photons, 20 ˝C, 1000 ˝C¨days, 100%, 20 ˝C, 1000 ˝C days, 0, 0, 0.5 for P opt , α opt , T opt , D opt , β h , β T , β D , γ h , γ T , and k rub , respectively.All 2094 half-hourly data points were used in optimization.A test with bootstrapping (20 iterations) to obtain parameter uncertainties for cases two and three resulted in no significant change in parameter estimates (results not shown), while bootstrapping could not be performed for case four due to hardware limitations.Figure 3 provides an overview of the modeling design and embedded test functionality to assess the verisimilitude of the model and effects of model approximations.
Remote Sens. 2015, 7, page-page 8 convergence, χ was constrained to 0.9, which is slightly higher than 0.8 used by [17] yet within the range specified by [59] for leaf-level photosynthesis (i.e., 0.5 to 0.9).In addition, some parameter values used in cases three and four were fixed to estimates obtained from cases two and three, respectively.Temperature parameters (i.e., , , ) were fixed in case three and both temperature and relative humidity parameters (i.e., , , , , ) in case four.Initial values for free parameters were 10 μmol•C•m −2 •s −1 , 0.1 μmol•C•μmol -1 photons, 20 °C, 1000 °C•days, 100%, 20 °C, 1000 °C days, 0, 0, 0.5 for , , , , , , , , , and , respectively.All 2094 half-hourly data points were used in optimization.A test with bootstrapping (20 iterations) to obtain parameter uncertainties for cases two and three resulted in no significant change in parameter estimates (results not shown), while bootstrapping could not be performed for case four due to hardware limitations.Figure 3 provides an overview of the modeling design and embedded test functionality to assess the verisimilitude of the model and effects of model approximations.Table 1.Overview of the different model inversion cases, the down-regulation mechanisms investigated and their fixed and free parameters.

Case
Enabled Modifier Functions Fixed Parameters Free Parameters 1 None ,

Results
The geometrically explicit model of canopy structure provided estimates of E al by numerically integrating down-welling diffuse and direct PAR. Figure 4 shows a 2D histogram of E al for elements p l within 5 m of the plot centre vertical axis, for a clear and a cloudy day in May 2009, using GF K = 15%.The graph shows a large variation in E al throughout the canopy and negligible absorption of radiation in the lower canopy for which two main causes can be identified: First, the forest plot was modelled using template crowns that exhibit no decay of foliage due to shading.However, when trees grow in stands, some leaves fall below the compensation point where gross CO 2 uptake equals respiration and foliage is shed; this phenomenon was not included in scene construction (Section 2.3).Second, a more technical concern, relates to the sampling density of photons in the Monte Carlo ray tracing simulations.Accurate sampling of light absorption in the lower canopy requires the simulation of a larger number of photons; however, as the larger contribution to P gs is made by the upper third of the canopy, the effects of restricted sample sizes on P gs estimates are limited [21,59].

Results
The geometrically explicit model of canopy structure provided estimates of by numerically integrating down-welling diffuse and direct PAR. Figure 4 shows a 2D histogram of for elements within 5 m of the plot centre vertical axis, for a clear and a cloudy day in May 2009, using = 15%.The graph shows a large variation in throughout the canopy and negligible absorption of radiation in the lower canopy for which two main causes can be identified: First, the forest plot was modelled using template crowns that exhibit no decay of foliage due to shading.However, when trees grow in stands, some leaves fall below the compensation point where gross CO2 uptake equals respiration and foliage is shed; this phenomenon was not included in scene construction (Section 2.3).Second, a more technical concern, relates to the sampling density of photons in the Monte Carlo ray tracing simulations.Accurate sampling of light absorption in the lower canopy requires the simulation of a larger number of photons; however, as the larger contribution to is made by the upper third of the canopy, the effects of restricted sample sizes on estimates are limited [21,59].Figure 5 shows the flux-tower derived , simulated , without the use of modifier functions for down-regulation (i.e., case 1) and using = 15%.The scatter plot shows a positive bias among a selection of the simulation results, and an overall low fraction of explained variation (r 2 = 0.66).To investigate this bias, data points were stratified by and ℎ (Figure 6).Plots along the horizontal are stratified by and along the vertical by ℎ.For and ℎ close to the optimum for photosynthesis, i.e., towards cooler temperatures and a relative humidity between 75% to 100%, a clear linear relationship between simulated and eddy-covariance derived was observed, whereas with distance from this meteorological optimum the relationship between simulated and measured productivity weakens and eventually becomes absent (i.e., for around 25 °C and ℎ < 50%).Figure 5 shows the flux-tower derived P gs,EC versus simulated P gs,sim without the use of modifier functions for down-regulation (i.e., case 1) and using GF K = 15%.The scatter plot shows a positive bias among a selection of the simulation results, and an overall low fraction of explained variation (r 2 = 0.66).To investigate this bias, data points were stratified by T and h (Figure 6).Plots along the horizontal are stratified by T and along the vertical by h.For T and h close to the optimum for photosynthesis, i.e., towards cooler temperatures and a relative humidity between 75% to 100%, a clear linear relationship between simulated and eddy-covariance derived P gs was observed, whereas with distance from this meteorological optimum the relationship between simulated and measured productivity weakens and eventually becomes absent (i.e., for T around 25 ˝C and h < 50%).  1 and text for an overview of parameters used for the different cases and Table 2 for an overview of parameter estimates obtained through model inversion.1 and text for an overview of parameters used for the different cases and Table 2 for an overview of parameter estimates obtained through model inversion.1 and text for an overview of parameters used for the different cases and Table 2 for an overview of parameter estimates obtained through model inversion.Figure 7 shows the behaviour of the modifier functions (y-axis) as they deviate from the optimum environmental condition (e.g., optimal temperature or relative humidity) expressed along the x-axis, and a scaling parameter that controls the pace of down-regulation with distance from the optimal condition.Using the modifier functions T mod and h mod (i.e., case 2), the photosynthetic capacity was found optimal at humid conditions (h > 75%) and a temperature of 19.5 ˝C.With the exception of χ, the duration of the initial linear portion of the light response curve, the optimization was unconstrained and photosynthetic capacity and quantum yield were obtained within their reported ranges [55,59].For example, the quantum yield α opt was found to be 0.097 (Table 2) which is among the lower range (i.e., 0.09 to 0.11 µmol¨C¨µmol ´1 photons) reported by [55] for the leaf level, while values for optimal photosynthetic capacity are commonly found between 11 and 27 µmol C m ´2¨s ´1 for the leaf level (e.g., Jones 1992 in [59]).Figure 8 demonstrates improvements in using the modifier functions (cases 2-4) for the same GF K used in Figures 4-6 (i.e., GF K = 15%).The scatter plot (Figure 8a) shows a much improved correspondence with P gs,EC ; however, the results indicate that biases remain that are associated with the diurnal response averaged over the study period (Figure 8g) and the seasonal course of P gs,EC (Figure 8j).Underestimates at noon may indicate that light inhibition of foliage respiration is underestimated in EC data.The decline in P gs,EC shows a more linear course from 13:00 PST onward than P gs,sim .This may indicate that down-regulation of photosynthesis is more linear with changes in relative humidity than what is currently modelled using the bell-shaped modifier functions.For example, Makela et al. [11] have used a linear function and a threshold to model temperature effects on the light use efficiency, and used an exponential relationship with vapour pressure deficit.Seasonal trends in P gs can be expected and result from different phenological stages of growth.To capture seasonal variation, the third case was investigated where T mod parameters were fixed to estimates obtained in case two and parameters for h mod and D mod were unconstrained (Figure 8b,e,h,k).The correspondence with P gs,EC further improved from r 2 = 0.71 to r 2 = 0.74 and the results showed no significant bias between months.Inclusion of shade acclimation of leaves in case four caused a minor improvement in the correlation with P gs,EC (Figure 8c) and an improved precision at Noon and during the afternoon can be observed (Figure 8i), where presumably the distribution of radiation across the canopy matches the nitrogen allocation profile.Finally, Table 2 shows the sensitivity of photosynthesis parameters to changes in gap fractions, GF K .
Remote Sens. 2015, 7, page-page Figure 7 shows the behaviour of the modifier functions (y-axis) as they deviate from the optimum environmental condition (e.g., optimal temperature or relative humidity) expressed along the x-axis, and a scaling parameter that controls the pace of down-regulation with distance from the optimal condition.Using the modifier functions and ℎ (i.e., case 2), the photosynthetic capacity was found optimal at humid conditions (ℎ > 75%) and a temperature of 19.5 °C.With the exception of , the duration of the initial linear portion of the light response curve, the optimization was unconstrained and photosynthetic capacity and quantum yield were obtained within their reported ranges [55,59].For example, the quantum yield was found to be 0.097 (Table 2) which is among the lower range (i.e., 0.09 to 0.11 μmol•C•μmol −1 photons) reported by [55] for the leaf level, while values for optimal photosynthetic capacity are commonly found between 11 and 27 μmol C m −2 •s −1 for the leaf level (e.g., Jones 1992 in [59]).Figure 8 demonstrates improvements in using the modifier functions (cases 2-4) for the same used in Figures 4 to 6 (i.e., = 15%).The scatter plot (Figure 8a) shows a much improved correspondence with , ; however, the results indicate that biases remain that are associated with the diurnal response averaged over the study period (Figure 8g) and the seasonal course of , (Figure 8j).Underestimates at noon may indicate that light inhibition of foliage respiration is underestimated in EC data.The decline in , shows a more linear course from 13:00 PST onward than , .This may indicate that down-regulation of photosynthesis is more linear with changes in relative humidity than what is currently modelled using the bell-shaped modifier functions.For example, Makela et al. [11] have used a linear function and a threshold to model temperature effects on the light use efficiency, and used an exponential relationship with vapour pressure deficit.Seasonal trends in can be expected and result from different phenological stages of growth.To capture seasonal variation, the third case was investigated where parameters were fixed to estimates obtained in case two and parameters for ℎ and  for the absorption probability that estimates the fraction of photon flux passing through the reference plane from a specific direction that is absorbed by a specific scene element.The model was parameterized using 3D structural data that was acquired with laser scanning technology and field observations of crown architecture, time series measurements of incident diffuse and direct PAR, and leaf-level photosynthesis parameters.While the virtual scene is constructed in forward modeling mode, the photosynthetic parameters were determined through model inversion against EC-derived P gs .Despite differences in scale between the modeled forest plot and the footprint of the EC flux tower, the model optimization resulted in simulated P gs values matching the range of error typical for the EC technique (e.g., [12,17,60]).Moreover, photosynthesis model parameters P opt and α opt were found within typical ranges for the leaf level, despite the coarse fidelity that was used to represent the canopy structure.When sun acclimation was introduced into the model using the nitrogen allocation parameter, k rub , estimates of P opt exceeded those reported in literature, suggesting that other acclimation models need to be considered (e.g., [61]).Our finding that diffuse E ic is a major control on P gs confirms the findings by [62,63].For example, Cai et al. [63] found that P gs of this stand can be effectively modelled with a single big leaf model using the sum of down-welling diffuse E ic and a fraction (of approximately 20%) of down-welling direct E ic .
A number of limitations of the proposed model remain, including neglecting the understory vegetation, the fidelity of the three-dimensional forest scene, assumptions of homogeneity of relative humidity and temperature within the canopy volume.The modelled scene does not represent individual needles or the clumping of needles into shoots, but rather comprises a limited number of canopy elements to represent foliage.Shading within shoots establishes significant effects on quantum yield and photosynthetic capacity [17,64] that are excluded using our coarse-resolution scene.Inclusion of finer levels of detail is no restriction to the ray tracing method, however, and is increasingly supported through advances in (mobile) laser scanning technology and other 3D acquisition techniques.Widlowski et al. [20] have demonstrated the use of ray tracing for scenes consisting of several millions of elements and for scenes comprising a stand of individual trees at needle-level resolution.For cases where the data size exceeds computer memory and the scenes exhibit some repetition in structure (e.g., scenes that represent fractal geometries or a collection of identical objects), techniques known as instantiation exist to reduce memory requirements (e.g., [35]).To ensure that absorption probabilities at the level of individual scene elements are accurately sampled, a higher structural fidelity of scenes requires the simulation of a larger numbers of photons and requires the use of advanced, optimized ray tracing algorithms.Increasing scene complexity in ray tracing simulations also provides for studies that investigate the physical mechanisms behind empirical observations of canopy reflectance and its linkages with leaf optical properties [65] and can further be used to analyze the distribution of radiation within shoots and needles [64].
While the current study focused on light transport, no links were established between shortand long-wave radiation and evapotranspiration.Forest canopies provide a level of insulation to extreme weather conditions and establish more temperate, humid conditions closer to ground level [66].The current model does not simulate heat transfer with canopy depth or the effects of latent heat transfer resulting from evapotranspiration, causing overestimates of drought and temperature stresses on the lower parts of the canopy.The effects on stand-level gross primary productivity may be relatively small, however, since the largest photosynthetic contributions come from the upper canopy (e.g., [59]).Efforts to include thermal energy transfer may incorporate radiosity modeling [25,52] or using 1D radiative transfer equations, acknowledging that the variability of temperature within the canopy space is relatively small [61] compared to light gradients [67].In addition, no leaf-level measurements of photosynthetic parameters were made in this study, and as a result, no validation or inference around the variability of photosynthetic parameters within the canopy space can be made at this stage.For comparison, Ibrom et al. [17] estimate that physiological variability among needles contributes to over 6% of measurements of electron transport rate-i.e., analogous to measurements of quantum yield-owing to differences caused by factors including shoot morphology and self-shading, leaf acclimation [61] and changes in hydraulic conductance with height [68].In addition, the modifier functions did not consider hysteresis to changes in temperature and relative humidity.It is likely, however, that photosynthesis responds more quickly to precipitation than to dry conditions.For example, for a drought to have a significant impact on photosynthesis, a weather record of several weeks may need to be considered, whereas a sudden event of heavy rain after a prolonged period of dry conditions affects the growth of vegetation more quickly.Other important considerations that are currently omitted in our model are the co-limitation effect of nitrogen availability on electron transport and photosynthetic capacity [61], as well as the acclimation of leaves to seasonal temperature fluctuations [69] and the effects of transient changes in incident radiation on photosynthetic down-regulation [70].
Whereas the parameter values retrieved in this study coincide with typical ranges reported in literature, caution should be taken towards inference from model inversion, since multiple combinations of parameter values can result in similar model estimates: an effect known as equifinality [17] and this particularly holds when larger numbers of parameters are estimated; however, the finding that most parameter estimates were obtained within previously reported ranges at the leaf level demonstrates a potential for the combined use of leaf and stand-level measurements to gain an understanding in the relative importance of physiological processes that establish the observed eddy-covariance fluxes [17].Whereas several one-dimensional models exist that produce similar accuracies against eddy-covariance flux measurements, these models are ultimately limited in the representation of processes at scales finer than the stand level.Three-dimensional plot-level models address the fusion of leaf and stand-level CO 2 exchange and consider a wide range of physiological processes.The finer spatial and temporal resolutions covered by these models enhance the opportunity to couple leaf and stand-level measurements such as obtained from cuvette and leaf chambers [17] and leaf spectral information [71].An effective way to combine shoot and canopy level information may be through the use of proximal narrow-waveband sensors data [72] installed within eddy-covariance footprints.These sensors provide year-round information about energy dissipation through the xanthophyll cycle: a predominant mechanisms through which reversible down-regulation of photosynthesis is achieved [73,74].Augmenting plot-level models using these sensor data may provide constraints on parameter ranges to avoid or reduce effects of equifinality, and can serve to validate the simulated radiation regime at specific locations within the canopy.Ultimately, the knowledge gained from plot-level P gs models and in situ sensor network data may be used to formulate empirical relationships that can be used to calibrate and validate one-dimensional P gs models that operate over broader scales.

Conclusions
Canopy structure is an important driver for canopy PAR modeling and is an important component in forest gross primary productivity (P gs ) models as it regulates how portions of diffuse and direct PAR are distributed throughout the canopy.In this study, a method was presented that provides for rapid updates of the canopy radiation regime, after an initial, computationally intense model initialisation phase.The model relies on the population of a look-up table that stores probabilities of absorbing PAR by individual scene elements and handles the dynamic changes in the canopy PAR that result from solar tracking and changes in atmospheric conditions.Once canopy radiation was simulated, a single leaf-level photosynthesis model was used to model canopy productivity for a time series of over 2000 measurements spread throughout a growing season.Through optimization of the model parameters to EC-derived P gs it was found that photosynthetic assimilation was highest at 19 ˝C and humid conditions and during the mid-season, in agreement with reported values for this climatic zone.Improvements to the model representation, canopy structure, and the inclusion of the transfer of long-wave radiation, sensible and latent heat, and in situ sensor network data have considerable potential to refine these estimates and may lead to a tighter coupling between eddy-covariance stand level estimates of gross primary productivity and leaf-level observations of photosynthesis.Such coupling can be used to gain insight in function-structure relationships at spatial and temporal scales where laboratory and field knowledge of photosynthesis meet.
Acknowledgments: Parts of this research were funded with the NSERC Discovery grant to Nicholas Coops.Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors an do not necessarily reflect the views of the funding body.We are thankful to Robert Guy, Thomas Hilker and Mike Wulder for their insightful comments and support throughout the research, Darius Culvenor and Glenn Newnham from the Commonwealth Scientific and Industrial Research Organisation (CSIRO) for their help and lending of equipment for data acquisition, members of the Biometeorology and Soil Physics Group, Faculty of Land and Food Systems (University of British Columbia) for maintenance of the research site.
Author Contributions: Van Leeuwen designed and developed the simulation code, analyzed the simulation results, and provided the initial draft of the manuscript; Coops was involved as Van Leeuwen's advisor; Dr. Black provided the EC data along with key inputs as to the interpretation and analysis of these data; All co-authors were involved with the various stages of manuscript revision.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.An illustration of the level of detail in Arbaro tree modelling software.

Figure 2 .
Figure 2. (A) Forest plot reconstructions representing individual stems; and (B) a coarser level of detail, where foliage clustered around branches is represented as planar polygons.

Figure 1 . 5 Figure 1 .
Figure 1.An illustration of the level of detail in Arbaro tree modelling software.

Figure 2 .
Figure 2. (A) Forest plot reconstructions representing individual stems; and (B) a coarser level of detail, where foliage clustered around branches is represented as planar polygons.

Figure 2 .
Figure 2. (A) Forest plot reconstructions representing individual stems; and (B) a coarser level of detail, where foliage clustered around branches is represented as planar polygons.

Figure 3 .
Figure 3. Schematic diagram of methods including test functions; GPP = gross primary productivity, PAR = photosynthetically active radiation.

Figure 3 .
Figure 3. Schematic diagram of methods including test functions; GPP = gross primary productivity, PAR = photosynthetically active radiation.

Figure 4 .
Figure 4. Simulation results of absorbed PAR for 6:00 PST (a) and 13:00 PST (b) on the cloudy day of 3 May 2009; and for 6:00 PST (c) and 16:00 PST (d) on a sunny day (8 May 2009) plotted as a 2D heat map.The gap fraction,, of photosynthetic scene elements, , was set to 15%.The x-axis is incident PAR (μmol photons m −2 •s −1 ) and the y-axis is height above the ground level at the plot centre.Colours represent observation frequency for an arbitrary, constant bin size.

Figure 4 .
Figure 4. Simulation results of absorbed PAR for 6:00 PST (a) and 13:00 PST (b) on the cloudy day of 3 May 2009; and for 6:00 PST (c) and 16:00 PST (d) on a sunny day (8 May 2009) plotted as a 2D heat map.The gap fraction, GF K , of photosynthetic scene elements, p l , was set to 15%.The x-axis is incident PAR (µmol photons m ´2¨s ´1) and the y-axis is height above the ground level at the plot centre.Colours represent observation frequency for an arbitrary, constant bin size.

Figure 5 .
Figure 5. Model inversion for case 1 (2094 observations, no cross validation): absence of downregulation functions.Shown is (a) the correlation with eddy-covariance derived gross primary productivity estimates; (b) the distribution of measured and simulated estimates of gross primary productivity against incident photosynthetically active radiation; (c) average diurnal responses.See Table1and text for an overview of parameters used for the different cases and Table2for an overview of parameter estimates obtained through model inversion.

Figure 6 .
Figure 6.A stratification of the data points shown in Figure 5 into temperature and relative humidity classes.

Figure 5 .
Figure 5. Model inversion for case 1 (2094 observations, no cross validation): absence of down-regulation functions.Shown is (a) the correlation with eddy-covariance derived gross primary productivity estimates; (b) the distribution of measured and simulated estimates of gross primary productivity against incident photosynthetically active radiation; (c) average diurnal responses.See Table1and text for an overview of parameters used for the different cases and Table2for an overview of parameter estimates obtained through model inversion.

Figure 5 .
Figure 5. Model inversion for case 1 (2094 observations, no cross validation): absence of downregulation functions.Shown is (a) the correlation with eddy-covariance derived gross primary productivity estimates; (b) the distribution of measured and simulated estimates of gross primary productivity against incident photosynthetically active radiation; (c) average diurnal responses.See Table1and text for an overview of parameters used for the different cases and Table2for an overview of parameter estimates obtained through model inversion.

Figure 6 .
Figure 6.A stratification of the data points shown in Figure 5 into temperature and relative humidity classes.

Figure 6 .
Figure 6.A stratification of the data points shown in Figure 5 into temperature and relative humidity classes.

Figure 7 .
Figure 7.A graph illustrating the behaviour of the modifier function.The modifier function is a bell shaped curve whose response is constrained to the range 0 to 1 (y-axis).

Figure 7 .
Figure 7.A graph illustrating the behaviour of the modifier function.The modifier function is a bell shaped curve whose response is constrained to the range 0 to 1 (y-axis).

Table 1 .
Overview of the different model inversion cases, the down-regulation mechanisms investigated and their fixed and free parameters.

Table 2 .
Overview of parameter estimates obtained through model inversion.
˝C¨days sensitivity of P max to changes in D β E µmol¨photons¨m ´2¨s ´1 sensitivity of P max to changes in E al β T ˝C sensitivity of P max to changes in T γ h lag in response of P max to changes in h γ T lag in response of P max to changes in T Where the dimensions are not specified, the quantities are dimensionless. *