Competition for light between individual trees lowers reference canopy stomatal conductance : Results from a model

[1] We have used an ecosystem model, TREES (Terrestrial Regional Ecosystem Exchange Simulator), to test the hypothesis that competition for light limits reference canopy stomatal conductance (GSref; conductance at 1 kPa vapor pressure deficit) for individual tree crowns. Sap flux (JS) data was collected at an aspen‐dominated unmanaged early successional site, and at a sugar maple dominated midsuccessional site managed for timber production. Using a Monte Carlo approach, JS scaled canopy transpiration (EC) estimates were used to parameterize two versions of the model for each tree individually; a control model treated trees as isolated individuals, and a modified version incorporated the shading effects of neighboring individuals on incident radiation. Agreement between simulated and observed EC was better for maple than for aspen using the control model. Accounting for canopy heterogeneity using a three‐dimensional canopy representation had minimal effects on estimates of GSref or model performance for individual maples. At the Aspen site the modified model resulted in improved EC estimates, particularly for trees with lower GSref and more shading by neighboring individuals. Our results imply a link between photosynthetic capacity, as mediated by competitive light environment, and GSref. We conclude that accounting for the effects of canopy heterogeneity on incident radiation improves modeled estimates of canopy carbon and water fluxes, especially for shade intolerant species. Furthermore our results imply a link between ecosystem structure and function that may be exploited to elucidate the impacts of forest structural heterogeneity on ecosystem fluxes of carbon and water via LiDAR remote sensing.


Introduction
[2] In forest ecosystems the ability to accurately quantify radiation budgets in heterogeneous canopies is a key challenge associated with modeling carbon and water fluxes [Muraoka and Koizumi, 2005].This is particularly important for land surface models (LSM) and dynamic global vegetation models (DVGM) that are used intensively to understand ecosystem responses and feedbacks to climatic change [Asner et al., 1998].Several solutions of varying complexity exist to deal with this problem.The simplest method for modeling canopy processes related to stomatal function is to treat the entire canopy as a single leaf.Although the simplicity of these big-leaf models is attractive, it can lead to insufficient model accuracy [Baldocchi et al., 1991;Raupach and Finnigan, 1988], particularly in heterogeneous ecosystems.A clumping parameter (W) can be employed to overcome the assumption of randomly distributed foliage associated with measures of leaf area index (L) incorporated in such models [Chen et al., 2008;Nilson, 1971].A second approach differentiates between sunlit and shaded canopy layers, effectively treating the canopy as two big leaves [Depury and Farquhar, 1997].Horizontal heterogeneity can be accounted for with W, and accounting for differences between sunlit and shaded canopy layers incorporates vertical heterogeneity.Although it fails to capture differences in incident radiation resultant from vertical heterogeneity among individuals, this multilayered approach is often seen as a positive tradeoff between accuracy and model complexity.A third more complex method involves representing individual crowns in a threedimensional array and calculating explicit ray traces through the canopy space between the sun and a point (crown) of interest [Brunner, 1998;Cescatti, 1997].Obviously this approach is limited by data requirements.However, it offers an explicit representation of canopy heterogeneity, and so is typically more accurate.Furthermore, three-dimensional models enable simulation studies in which canopy architecture can be manipulated to understand its effects on radiative transfer, as well as studies that explore succession and ecosystem demography [Courbaud et al., 2003;Genard et al., 2000].As such these models are relatively popular, despite their high complexity and computational requirements, especially for agro-forestry applications [Bartelink, 1998].
[3] Plasticity in photosynthetic capacity within species, expressed via changes in leaf structure [Sack et al., 2006], leaf mass [Sims and Pearcy, 1994], or biochemical limitations [Lei and Lechowicz, 1997a;Naidu and Delucia, 1997], as a response to changes in radiation environment is well documented.Additionally, there is mounting evidence that photosynthetic capacity and plant hydraulics are linked for a wide range of species and plant functional types [Allen and Pearcy, 2000;Brodribb and Feild, 2000;Brodribb et al., 2002;Campanello et al., 2008].From a modeling perspective these factors are most important when considering individual crowns, and become less important with increasing scale, as ecosystem fluxes become scale invariant [Enquist et al., 2007].However, functional differences expressed at the canopy scale associated with variation in intercepted radiation are important for parameter estimation and moving across scales [Mackay et al., 2010].
[4] At two forest stands in northern Wisconsin (shadetolerant sugar maple (Acer saccaharum Marsh) and shadeintolerant trembling aspen (Populus tremuloides Michx) dominated), spatial variability in transpiration on an individual tree basis has been characterized and partially explained by tree size [Loranty et al., 2008;Traver et al., 2010].Further, we have shown that ignoring local variability in stand structure can lead to errors in estimates of stand transpiration scaled from sap flux observations [Mackay et al., 2010].An inverse relationship between stomatal conductance and a canopy competition index (C I ) further explains differences in stomatal control of water loss between individuals [Loranty et al., 2010].The following equation was used to calculate C I : where for each competitor i, a is the degree of southness to the sap flux tree when a ≤ 180°and ½ azimuth direction where a > 180°, d i is the distance between the sap flux tree and the competitor, h s is the height of the sap flux tree, and h i is the height of the competitor.This metric is unique from other weighted distance measures of competition because it includes information on vertical heterogeneity via the inclusion of height.Independent analyses of height, distance, and direction revealed little or no correlation with stomatal conductance indicating that each component of C I is important and suggesting that the degree of competitive shading was responsible for variability in stomatal conductance between individual trees.The observed relationship was stronger for aspen than it was for maple.Thus, we hypothesize that this relationship is a consequence of competition for light between individual tree crowns, and that this competition contributes to variability in stomatal conductance through light-limited photosynthesis.Additionally we hypothesize that the magnitude of variability in stomatal conductance attributed to competitive shading will be greater in aspen than maple because maple is a more shade tolerant species.In the current study we test these hypotheses by incorporating a three-dimensional canopy representation into an ecosystem model that constrains stomatal conductance using a combination of plant hydraulics and photosynthesis.

Study Site
[5] The study was conducted near the town of Park Falls (45.95 N, 90.27 W), in northern Wisconsin in the Chequamegon-Nicolet National Forest, at two sites located less than 2 km from the WLEF eddy flux tower [Bakwin et al., 1998] as part of the Chequamegon Ecosystem Atmosphere Study (ChEAS).The area represents the interface between northern temperate and southern boreal ecosystems and is part of the northern highlands physiographic province.Geomorphic features in the area are outwash, pitted outwash, and moraines resulting in gently rolling topography.Climate is characterized by long winters and a short growing season with mean January and July temperatures of −12°C and 19°C respectively [Fassnacht and Gower, 1997].
[6] Field observations for the study were collected at two sites during the summer of 2005.One site, hereafter referred to as Aspen, was a transition from forested wetland to upland, with the wetland dominated by speckled alder (Alnus incana Rugoi) and white cedar (Thuja occidentalis), and upland dominated by trembling aspen [Loranty et al., 2008].The second site was a sugar maple dominated stand that contained several patches of plantation red pine (Pinus resinosa), and is hereafter referred to as Maple [Traver et al., 2010].The stand at the Aspen site was regenerating from a clear-cut that occurred approximately 20 years ago, and subject to no further management.Little local relief existed at the Aspen site, with a maximum elevation difference of 3 m across the site resulting from a gentle upward slope from the wetland in the northwest moving to the south and east.Trees at the Maple site ranged from approximately 50 to 75 years in age.The site was previously managed for red pine timber production, and the maple portions are recruits that have established subsequently.The presence of stumps throughout the site suggests it has been thinned in the past and during the study there were trees flagged for removal.This suggests that management is on going although the exact prescription is not known.Elevation varied by approximately 15 m across the Maple site as a result of local microtopography.At each site a series of plots, with diameters of 5 m and 6.5 m at the Aspen and Maple sites respectively, were established using a cyclic sampling design [Burrows et al., 2002].Within each plot the largest tree from the dominant species was selected for sap flux instrumentation, with additional trees sampled where resources allowed.At each site tree size varied spatially [Loranty et al., 2008;Traver et al., 2010], so our sampling strategy effectively captured the entire range of tree size and competition at each site [Loranty et al., 2010].Approximately 200 individual aspen and maple were instrumented in total; however because variation in transpiration has been explained by tree size [Loranty et al., 2008;Traver et al., 2010] and because data were limited by the subset of days we chose to use for analysis (described below) 20 trees from each site were selected for analysis in the present study.

Sap Flux and Environmental Measurements
[7] Trees were instrumented for sap flux using the heat dissipation method [Granier, 1987].At the Aspen site measurements were taken from 27 May to 8 July and from 19 July through 19 August at the Maple site.Diameter at breast height (DBH) and local xy coordinates were recorded for each tree selected for sap flux instrumentation at both sites.Temperature and relative humidity (Vaisala HMP 45C, Vaisala Oyj, Helsinki, Finland) were measured at two-thirds canopy height at both sites.Agreement between in situ vapor pressure deficit (D) and 30 m observations from the WLEF tower indicate that the canopy was well coupled with the atmosphere [Ewers et al., 2008].Sap flux, temperature, and relative humidity measurements were recorded every 30 s (CR10X, Campbell Scientific, Logan, UT, USA) and aggregated to 30 min values.Wind speed, photosynthetically active radiation (Q 0 ), and precipitation measurements from the nearby WLEF (∼2 km) [Davis et al., 2003] and Lost Creek (∼10 km) [Cook et al., 2004] flux towers were used as model inputs.Measurements of leaf gas exchange were made at the Aspen site during the study period for a subset of 12 aspen trees using an LI 6400 (Li-COR Biosciences, Lincoln, NE).A scaffolding was used to sample six trees in the upland and six in the wetland portion of the site, and photosynthetic light response curves were generated using standard techniques [Long and Bernacchi, 2003].We were unable to access the canopy at the Maple site and used light response curves from a maple stand in the region (Eric Kruger, unpublished data).
[8] Transpiration per unit ground area (E C ) was estimated from sap flux measurements using the following equation: where J S is sap flux per unit xylem area, A S is sapwood area, A G and is ground area.A S was calculated from DBH using species specific allometric relationships established at nearby sites [Ewers et al., 2002].This previous study examined [Ewers et al., 2002] radial and circumferential trends in sap flux for each species in addition to measuring sapwood depth, and so we are confident that any observed differences between species are not a result of error associated with scaling from sensor to whole tree.A G was assumed to be the projected cross-sectional crown area for each individual [Kostner et al., 1992].Measurements taken as part of a stand inventory to characterize the canopy light environment (described below) were used to develop an allometric relationship between DBH and A G at the Aspen site.To accomplish this, two measurements of canopy diameter were used to calculate cross-sectional crown area, which was regressed against DBH and fit with a logarithmic equation.At the Maple site a similar approach was used, but the allometry was taken from the literature [Frelich and Lorimer, 1985].

Stand Measurements
[9] Stand inventories for each site were collected so that the canopy light environment could be modeled (Table 1).During August 2006 and August 2007, for the Aspen and Maple sites respectively, an inventory was recorded of all neighbors within 5 m and 10 m of each P. tremuloides and A. saccharum individual instrumented for sap flux, respectively.Trees lying to the north of the sap flux tree between 345°a nd 45°were excluded from the analysis on the grounds that their competitive effects within the canopy would be negligible since their crown was never collinear with the sun and crown of the sap flux tree.For each sampled tree local xy coordinates and DBH were recorded.At the Aspen site tree height (H T ) was measured for each tree using a laser rangefinder and clinometer by triangulation and for a subset of aspen trees crown diameter in two directions and height to base of live crown were measured so that allometric relationships could be established.At the Maple site a higher and denser canopy, a thick understory, and increased microtopography made height measurements impossible.Consequently tree heights at the maple stand were calculated using species-specific allometric equations developed for stands with similar structure and age in the region [Frelich and Lorimer, 1985;Perala and Alban, 1994;Ter-Mikaelian and Korzukhin, 1997].Differences in relative tree height at the Maple site were corrected using a LiDAR derived digital elevation model with 1 m resolution.The small elevation change across the Aspen site was monotonic and thus could be accounted for in the canopy model so no corrections were made to measured tree heights.
[10] The canopy light model (described below) required further biometric inputs of leaf area density (LAD), height to base of live crown (H C ), and canopy diameter, which were derived primarily from allometry.Table 1 shows the species makeup of each stand inventory.For the areas inventoried the species sampled for sap flux are clearly dominant, however other species were present.For some of these less common species allometry was scarce.Available allometry was closely scrutinized, as relationships are often very site specific [Gilmore, 2001].Subsequently, we aggregated several species on the grounds that allometry for the dominant species was more appropriate.Specifically, relationships developed for canopy diameter and H C for aspen at the Aspen site were applied to all deciduous species at the site.Species-specific allometry for maple, cherry, and birch in the region were developed at mature stands and were not suitable for the younger individuals at the aspen stand.Of the sap flux trees from the Aspen site used in this study only four had deciduous species other than aspen in the neighborhood sampled for competition, and of these five were cherry and three were birch, and only one of these was taller than the sap flux tree.At the Maple site allometry was available for maple, basswood, and red pine and these species accounted for over 95% of the sample populations.
For the remaining species allometry for sugar maple was applied.Although combined species allometry is not ideal the primary concern here is presence or absence of a neighboring crown competing for light.Moreover, any errors are likely negligible because of the relatively coarse resolution (2 m) of the three-dimensional canopy model that was employed to calculate path lengths through the canopy.

Experimental Approach
[11] A model-data fusion approach was used to examine the effects of canopy heterogeneity on canopy stomatal conductance at D = 1 kPa (G Sref ) for individual trees, ).E C from each of the resulting 10,000 simulations is compared to observed E C in order to identify the best simulation to use for analysis.In the light model (b) inputs from stand inventories are used to calculate the distance a ray of sunlight must travel through the neighboring canopy to reach the crown of interest using the tRAYci model.The TREES model is modified to accept path length as an input, which is used to constrain incident radiation and partition LAI into sunlit and shaded portions.All other inputs and outputs for TREES are the same, and an identical model evaluation procedure are followed to identify the best out of 10,000 simulations to be used for analysis.Note that unique estimates of LAI for each tree, based on allometry, are used as model inputs.
derived by inverting a canopy transpiration model (Figure 1).A control model (hereafter big leaf model) used a multilayer big-leaf canopy representation, and a modified model (hereafter light model) constrained the canopy light environment using a three-dimensional canopy model.A Monte Carlo technique was used for model parameterization, and an initial uncertainty analysis was used to constrain parameter ranges.Estimates of G Sref from the big leaf and light models were compared to assess the impacts of accounting for canopy heterogeneity.We chose to use data from cloudless days in order to more explicitly elucidate the effects of competitive shading.Cloudless days are defined as having peak Q 0 > 1800 mmol m −2 sec −1 and diel plots of Q 0 that exhibited minimal high-frequency microvariation in the form of spikes or other anomalies (indicative of intermittent cloud cover) (Figure 2).Data collected during the days of 21-25 June and 1-2 July at the Aspen site, and 21, 27, 30 July and 4-7 August at the Maple site was chosen for analysis (Figure 2).The limited number of cloudless days reduced data availability and so we selected representative subsample of approximately 20 individuals from each site based on data availability associated with random power supply and site maintenance constraints.Data from each tree were used to parameterize the big leaf and light models.

Transpiration Model
[12] Terrestrial Regional Ecosystem Exchange Simulator (TREES) [Ewers et al., 2008;Mackay et al., 2003;Samanta et al., 2007;Samanta et al., 2008;Samanta and Mackay, 2003] is a modified big-leaf canopy model that simulates transpiration using a coupled hydraulic and biochemical model of stomatal conductance.Allometry is used to generate a unique value of LAI for each tree, and this is partitioned into sunlit and shaded elements based upon light attenuation and a species-specific clumping parameter.An initial estimate of hydraulically limited stomatal conductance is constrained using a photosynthesis model [Depury and Farquhar, 1997] to yield a final value of stomatal conductance.This process is carried out in parallel for sunlit and shaded portions of the canopy.The following paragraphs first provide a general overview of the model structure, and a more detailed description of portions that were altered for this study to differentiate the big leaf and light models is included in Appendix A. Model evaluation and parameterization for this study is described in the subsequent section.
[13] The following equation is used to generate an initial value of hydraulic limited stomatal conductance [Oren et al., 1999]: where G Sref0 is reference canopy average stomatal conductance at D = 1 kPa, and m is a parameter that represents the sensitivity of G S to D. Among isohydric species that maintain a minimum leaf water potential despite declines in soil water potential [Tardieu and Simonneau, 1998] a 0.6 proportionality that has been theoretically explained in terms of plant hydraulics [Oren et al., 1999] has been observed across a broad range of species and conditions [Addington et al., 2004;Ewers et al., 2005;Ewers et al., 2007].
[14] Canopy element (k) leaf specific CO 2 conductances (g c0,k ) were calculated using the following equation: where L k is element leaf area index, g va is aerodynamic conductance taking into account turbulent and laminar heat and vapor conductance [Campbell and Norman, 1998], and 1.6 is the proportionality between molar H 2 O and molar CO 2 .Element hydraulically constrained stomatal conductances to CO 2 were subsequently used to calculate photosynthesis using a hybrid model [Depury and Farquhar, 1997;Katul et al., 2003]  ) is calculated as: where A n is net photosynthesis, c a is atmospheric CO 2 concentration, and c i,k is intercellular CO 2 concentration.This relationship assumes atmospheric CO 2 and leaf surface CO 2 concentrations to be equivalent.At this point canopy element transpiration was calculated using the Penman-Monteith Combination Equation [Monteith, 1965].Further description of TREES, including incorporation of the light model described below, can be found in Appendix A.

Light Model
[15] A spatially explicit three-dimensional canopy light model, tRAYci [Brunner, 1998], was used to estimate incident radiation for individual trees.tRAYci uses exact locations for all trees in a plot and geometric crown representation to model entire stands.A ray-tracing algorithm is used to calculate the path length that light rays travel through the canopy between the top of the canopy and a specified point of interest.From this calculation, the percent of above canopy light (PACL) that reaches the points can be calculated using a simple Beer's law extinction coefficient.Path length calculations for individual trees that included only neighboring trees (or excluded self-shading) were utilized for the present study.
[16] Model inputs required for tRAYci were xyz coordinates, H C , up to four crown radii, and species.Species-specific parameters were used to further constrain light transmission through the canopy.Because our primary concern is path length, only the parameters used to characterize species crown geometry in terms of shape and foliage distribution are described here.For each stand separate crown shapes were input for deciduous and coniferous species.The model uses a three-parameter ellipse following Koop [1989], that requires two parameters that represent the exponent of a quarter ellipse.These parameters dictate the shape of the top and bottom of the crown, and allow flat, concave, convex, and cone shaped surfaces to be represented.The third parameter specifies the height within the live crown at which maximum crown width occurs in percent.Using these parameters conifers at both sites were represented as cones, and deciduous species were represented as ellipsoids.Mean values of LAD calculated with allometry were used for each species type.The model was set to the highest ray resolution, so that a sample ray was traced for every degree of zenith and azimuth yielding a total of 32,400 sample rays.For each tree we used the sun track from the median day of the sample period for each site to determine path length for each 30 min time step.Path length estimates were used to constrain incident radiation required to calculate photosynthesis, which was used to constrain stomatal conductance in a transpiration model described in the previous section.A more complete description of equations used to calculate incident radiation for sunlit and shaded canopy portions, in addition to the modifications made to account for canopy heterogeneity can be found in Appendix A.

Model Parameterization and Evaluation
[17] The results from the light model are compared to results from the big leaf model, and methods used to parameterize each model and define G Sref were identical.Briefly, an initial set of simulations where model parameterizations of G Sref0 and m were evaluated by comparing simulated G S against an observed value of G S derived empirically by inverting the Penman-Monteith equation [Monteith and Unsworth, 1990] in order to define constrained parameter ranges using the Monte Carlo sampling approach described by Samanta and Mackay [2003].Next two sets of 10,000 simulations, using the big leaf and light models, were generated for each tree using the constrained parameter ranges for G Sref0 and m (Figure 1).Simulated values of E C were compared to sap flux scaled E C to evaluate model performance.Index of Agreement (IOA) is a robust metric of model performance useful for making comparisons between models [Willmott, 1982].We used IOA to identify the top simulation from each model for each tree, and IOA in addition to r 2 to compare models.For the best model for each individual we fit a line to the upper envelope of modeled G S and m values to derive estimates of G Sref from both models Additionally we examined the relationship between G Sref and m as check because prior studies have found a 0.6 proportionality (3) between these two variables at similar sites in the region [Ewers et al., 2007;Mackay et al., 2003].

Results
[18] Profiles of Q 0 on days selected for analysis indicate absence of cloud cover (Figures 2b and 2d).E C observations at the Aspen (Figure 2a) and Maple (Figure 2c) sites corre-spond with diurnal patterns of D and Q 0 , and are representative of the ranges of E C for each species.The range and absolute magnitude of fluxes observed among aspen was greater than those observed among individual maple.A higher degree of short-term variability was observed within diurnal fluxes at the Aspen site.
[19] Values of E C simulated for all trees with the big leaf and light models showed acceptable agreement with observations for both sites (Figure 3).Overall, agreement between observed and modeled E C was better at the Maple site the Aspen site.In comparison to the big leaf model the light model resulted in improved agreement between simulated and observed E C at the Aspen site (Figure 3) and a decrease in model performance at the Maple site (Figure 3).Model parameterization remained consistent between models in terms of convergence for IOA, G Sref0 , and m (Table 2).
[20] Incorporating canopy heterogeneity in TREES had varying effects on model performance within species.At the Maple site agreement between simulated and observed E C with the Big Leaf Model was high (IOA > 0.9) (Figure 3 and Table 2), particularly in comparison to aspen.Simulated values matched observations especially well early and late in the day, times when D is typically low and Q 0 is limiting, and deviations typically occurred during periods of peak fluxes near midday (Figures 4a, 4c, and 4e).Subsequently incorporating canopy heterogeneity with the light model resulted in a slight decrease in model performance.Compared to big leaf simulated E C , mean IOA decreased from 0.95 to 0.94 (Table 2) with the light model.Mean r 2 for individual trees decreased from 0.83 to 0.81 with the light model, and the magnitude of decrease was inversely related to G Sref , with slight improvements for some individuals with low conductance (Figure 4c).Reductions in A N with the  4b, 4d, and 4f).Corresponding reductions in modeled G S effectively reduced E C during midday periods of peak fluxes, in most cases (Figures 4a, 4c, and 4e), resulting in the overall decreases in model performance mentioned earlier.
[21] At the Aspen site the light model showed better agreement between observed and simulated E C in comparison to the big leaf model.Observed E C exhibited a higher degree of short-term variability than maple at diurnal timescales (Figure 2), and subsequently there were differences in both the magnitude and timing of diurnal E C fluxes simulated with the big leaf model (Figures 5a, 5c, and 5e).The mean IOA for the top model from each tree improved from 0.86 with the big leaf model to 0.90 with the light model, and this corresponded to an increase in mean r 2 from 0.59 to 0.69 (data not shown).Among individuals the magnitude of increase in model performance was not consistent, rather there was a nonlinear inverse relationship between G Sref and the change in r 2 (Figure 6c).For individuals with high G Sref values (Table 2; i.e., D7 and E7), r 2 was relatively similar between models.Conversely individuals with low G Sref , greater short-term temporal variability in E C , and a greater degree of competitive shading, exhibited greater increases in r 2 .Trees with high G Sref and little competition for light exhibited little or no reduction in A N for the sunlit portion of the canopy (Figure 5f), while trees with lower G Sref and greater competition for light showed decreases in A N for sunlit leaves (Figures 5b and  5d).This finding generally resulted in changes in both the timing and magnitude of E C that were in better agreement with observations (Figures 5a, 5c, and 5e), resulting in the previously described increases in model performance.
[22] Values of G Sref for individuals generated using the big leaf and light models exhibited relatively little variation between the two models in comparison to within species variation (Figure 6a).Increases in G Sref with the light model were observed primarily at the Aspen site (Figure 6a).The models also adequately captured the 0.6 proportional relationship between m and G Sref (3) (Figure 6b), and this is    and (d), and D7 (e) and (f).consistent with observations [Ewers et al., 2007] and simulations for these species [Mackay et al., 2003].
[23] Results from individual trees can be used to aid interpretation of model results.For example, note that for tree A11 at the Aspen site (Figure 5c) simulated E C matched the observations reasonably well and that the differences between E C for the different models were relatively small.The light response curve (Figure 5d) shows that A n is reduced, however not to the linear portion of the plot at low Q 0 that corresponds to light-limited photosynthesis.This indicates competitive shading has not reduced irradiance sufficiently to limit photosynthesis.At the Maple site, where simulated E C matched observed values well in comparison to Aspen (Figures 4, 5, and 6), key differences were less micro-variability in peak fluxes and a greater reduction in A n (Figure 2).Using tree E6 as an example it can be seen that simulated E C from the light model was lower than that of the big leaf model.There was a substantial reduction in A n with Q 0 (Figure 5d).However, there were no differences between G Sref for the big leaf and light models despite reduced A n .One possible explanation for this lies in the canopy structure of maple.Values of L are higher at the Maple site and the species-specific value of W = 1.0 is higher than W = 0.6 for the Aspen.According to (A1) and (A2) then, the proportion of L* (sunlit leaf area) at the Maple site was lower than at the Aspen site.The light model modifies only direct beam radiation (Q b ) and so simulated values of shade leaf A n simulated with big leaf and light models were relatively invariant.As such reductions in sunlit leaf A n have less impact on G Sref , especially during periods of high Q 0 .

Discussion
[24] This study incorporated canopy heterogeneity into an established ecosystem model to test the hypothesis that light-limited photosynthesis as a result of canopy heterogeneity contributes to variability in stomatal conductance between trees for two forest stands dominated by different species.The light model improved simulated E C for aspen, but not for maple (Table 2).Overall, our results support the hypothesis that competition for light between crowns as a result of canopy heterogeneity contributes to variability in stomatal conductance, and the assertion that the magnitude of this effect is smaller in species that are more shade tolerant.In the following sections we interpret model results in this context, and discuss the implications of our findings at a range of scales.

Interpretation of Model Results
[25] Our modeling approach relies on observed E C so accurate scaling of J S is crucial, as any scaling errors would influence the results.As previously described, the method used estimate A S and scale point measurements of J S to the whole tree were developed in a previous study at a nearby Figure 6.Reference stomatal conductance plotted against the difference in r 2 between the big leaf model and the light model for (a) aspen, and (b) maple.The dotted and solid lines represent significant relationships between variables on cloudless days for aspen with the light model (r 2 = 0.71) and big leaf (r 2 = 0.62) models respectively.site [Ewers et al., 2002].These methods account for speciesspecific radial and circumferential trends in J S related to xylem anatomy and sapwood depth.Additionally, simulated values of G Sref observed in this study are similar to previous observations for both species, [Ewers et al., 2005;Ewers et al., 2007].Moreover, variability in modeled G Sref supports the idea that increasing spatial variability in E C with increasing D, observed at these sites and others, [Adelman et al., 2008;Loranty et al., 2008;Traver et al., 2010] is explained by differences in G Sref between individuals.Aspen exhibited greater variability in G Sref than maple, which can be expected in this context given that the magnitude of spatial variability at high D was greater at the Aspen site than at the Maple site [Traver et al., 2010].Thus we are confident that differences in E C and G Sref between individuals and species, and consequently our results, are not due to errors associated with scaling.
[26] Mean r 2 for maple simulations with the light model decreased slightly in comparison to the big leaf model.However agreement between observed and simulated E C was quite high to begin with (Table 1), and despite the decrease in model performance maple light model simulations showed higher agreement with observations than did aspen light model simulations, on average (Table 1).As a shade tolerant species [Ellsworth and Reich, 1993] mature maple stands are typically characterized by dense, closed canopies, with relatively low light transmission to the understory [Tobin and Reich, 2009] that implies a high degree of shelf-shading.This suggests that mainly horizontal (two-dimensional) heterogeneity affects G Sref at the Maple site.
[27] In contrast, light model simulations resulted in improved agreement between simulated and observed E C at the Aspen site as indicated by IOA (Table 2) and r 2 (Figure 6c).Increases in model performance, bolstered by reductions in simulated A N (Figure 5b, 5d, and 5f), indicate that competitive shading caused by vertical canopy heterogeneity is a source for variability in canopy stomatal conductance for aspen.This assertion is further supported by the inverse relationship between the magnitude of increase in r 2 and G Sref (Figure 6c), the link between G Sref and its competitive environment [Loranty et al., 2010], and previous studies illustrating a within species coordinated adaptive response to light environment between photosynthetic capacity and stomatal conductance [Lei and Lechowicz, 1997b;Tinoco-Ojanguren and Pearcy, 1993].This overall improvement indicates that accurate characterization of the canopy light environment is an important factor for modeling stomatal conductance in aspen, which is an early successional shade-intolerant species [Roden and Pearcy, 1993] characterized by more heterogeneous open canopies with higher understory light transmission [Chen et al., 1997].
[28] The relative complexity of our modeling approach warrants consideration of additional sources of error.Tree crown data used to model the light environment at the Maple site were derived using allometry from the literature that was corrected for elevation differences using a high resolution DEM.Consequently the potential error associated with these estimates [Fehrmann and Kleinn, 2006;Keller et al., 2001] is undoubtedly higher than that for the Aspen site where the same data were based on site-specific observations with no need for elevation correction.Subse-quently overestimation of competitive shading by the light model, compounded by species differences in shade tolerance and growth strategy [Kubiske et al., 2002;Niinemets and Tenhunen, 1997], may partially explain the modest decrease in model performance.Generally lower rates of E C may also contribute to results for maple, as it is known that E C is relatively insensitive to changes in G S when E C is low [McNaughton and Jarvis, 1991].It follows then that changes in G S at the Maple site may be too small to discern because E C is low.Additionally maple has higher L, with a relatively high degree of self-shading, such that changes in the proportion of sunlit and shaded L, and subsequently absorbed Q 0 , associated with competitive shading were relatively small in comparison to aspen.It may then be the case that a combination of uncertainties associated with methods of observation and data availability is responsible for the decrease in model performance with the increase in model complexity at the Maple sites.
[29] The assertion that a dynamic response to irradiance can be resolved at half-hourly time scales, and that there would be no residual effects of light limited photosynthesis in the morning on stomatal conductance in the afternoon, is supported by observations of rapid stomatal response to Q 0 and D for a number of species at a variety of scales [Zweifel et al., 2002].Therefore, the boundary line of the relationship between G S and D, used to calculate G Sref for both big leaf and light models, would reflect primarily hydraulic limitation at high D despite photosynthetic limitation at lower D during certain time periods.Thus, differences in G Sref between big leaf and light models can be interpreted as adjustments made by the big leaf model in an attempt to account for light-limited photosynthesis.However this could also apply to increases in G Sref by the light model in certain instances, because it is possible that the amount of competitive shading is overestimated and so higher G Sref would be necessary to compensate (i.e., aspen tree D7 in Table 2) (Figures 2 and 5).

Implications for Physiology and Modeling
[30] Our results point toward a positive relationship between canopy photosynthesis and G Sref , for aspen, more so than maple.A number of recent studies have reported coordination between plant hydraulics and photosynthesis [Brodribb and Jordan, 2008;Brodribb and Feild, 2000;Campanello et al., 2008;Maherali et al., 2008].This coordination is typically characterized by a reduction in hydraulic capacity that corresponds to a reduction in photosynthesis, and has been observed for a wide range of species [Brodribb and Feild, 2000].Manipulation experiments indicate that this relationship is dynamic and highly plastic [Campanello et al., 2008], and offer evidence for the genetic variability required for such plasticity [Callaway et al., 2003;Maherali et al., 2008;Weinig, 2000aWeinig, , 2000bWeinig, , 2000c]].Additional studies have observed direct relationships between stomatal density, leaf thickness [Niinemets, 2007], chlorophyll/chloroplast content [Boonman et al., 2009], and irradiance during leaf development within species.Recently Brodribb et al. [2007] found an inverse relationship between leaf anatomy and leaf hydraulic conductance, and a positive relationship between leaf hydraulic conductance and the maximum rate of photosynthesis.Although we are not able to identify the specific mechanism, our results support a relationship between photosynthetic capacity and plant hydraulics at the canopy level.
[31] This link between G Sref and photosynthetic capacity, and the apparent dynamic link between D and ecosystem fluxes [Adelman et al., 2008;Traver et al., 2010], modulated by differences in G Sref , could be particularly useful for quantifying flux responses to climate mediated changes in environmental and biological drivers [Bond-Lamberty et al., 2009].Indeed, it is hypothesized that declines in boreal productivity [Goetz et al., 2005] can be attributed to stomatal down regulation in response to increased D [Zhang et al., 2008].Our results along with those of other recent studies [Novick et al., 2009] suggest that variability in G Sref is linked to ecosystem structure, and thus, may potentially be observed with remote sensing tools such as LiDAR.However, the mechanistic explanation for this relationship is unclear as yet, and it is necessary to understand the processes driving ecosystem structure before ecosystem structure and function can be linked temporally.
[32] Despite their relative costs in terms of data requirements, three-dimensional canopy representations are necessary for understanding potential links between ecosystem structure and function [Genard et al., 2000;Lappi and Stenberg, 1998;Roupsard et al., 2008].Our results imply an apparent mechanistic link between ecosystem structure, photosynthetic capacity, and stomatal control of water loss for individual crowns.Understanding relationships between ecosystem structure and function is essential for more accurate scaling from points to areas [Mackay et al., 2010].Establishing such links has been recognized as an essential component for advancing regional and continental scale modeling studies of the biosphere, particularly in the face of environmental change [Medvigy et al., 2009;Moorcroft, 2003].High-resolution remote sensing platforms such as LiDAR capable of characterizing vegetation structure are increasingly prevalent [Brandtberg et al., 2003;Houldcroft et al., 2005], and efforts have been made to initialize ecosystem models with remote sensing based indices of vegetation structure Hurtt et al., [2004].Our results indicate that even multilayer big-leaf models may be insufficient for representing canopy fluxes at large spatial scales and that the amount of structural information required to sufficiently model canopy processes increases with the shadeintolerance of dominant species.

Conclusions
[33] Accounting for canopy heterogeneity by incorporating a 3-D canopy light transmission in a transpiration model improved model performance most in a shade intolerant species (aspen).For individual aspen growing in northern Wisconsin, the light model quantifies a link between plant hydraulics and light limited photosynthesis.In a shade tolerant maple stand, canopy competition for light appears to be less important as model performance with a multilayer big-leaf model was high, and inclusion of 3-D light transmission offered no improvement.However, this does not discount a link between plant hydraulics and photosynthetic capacity.Our results corroborate a link between within species difference in photosynthetic capacity imposed by ecosystem structure, and stomatal function.Future work should seek to further understand these links with the ulti-mate goal of improving large-scale estimates of ecosystem fluxes through platforms such as LiDAR remote sensing while simultaneously minimizing the amount intensive field observations required to accurately parameterize large scale models of ecosystem carbon and water fluxes.

Appendix A
[34] For light limited photosynthesis, absorbed radiation for sunlit and shaded canopy elements was calculated as follows.Sunlit leaf area of the canopy (L*), as described by Campbell and Norman [1998], is defined as: where L is leaf area index, and K be (Y) is the extinction coefficient for light in the canopy at zenith angle Y, based upon an ellipsoidal leaf angle distribution.The shaded portion of the canopy is then calculated by subtracting L* from L. K be (Y) is calculated as follows: where x is the leaf angle distribution, and W(y) is a canopy clumping parameter ranging from 0-1 that varies with zenith angle according to the following equation: where W is the clumping parameter when the canopy is viewed from nadir, and p is a variable related to the ratio of crown height to crown diameter.Typically x, W, and p are species-specific model input parameters derived from literature values [Campbell and Norman, 1998].
[35] Incoming solar radiation above the canopy, a model input, is partitioned into incoming direct (Q ob ) and incoming diffuse (Q od ) radiation as described in Spitters [1986].Within canopy total beam (Q bt ), direct beam (Q b ), and diffuse (Q d ) radiation fluxes are calculated using the following equations: where Y is the zenith angle of the sun, and t bt , t b , and t d is the proportion of total beam, direct beam, and diffuse radiation transmitted through the canopy.The following equations are used to calculate t bt , t b , and t d , respectively: where a is the absorptivity of leaves for radiation.Each component of incoming radiation described above is computed in parallel for photosynthetically active radiation (PAR), and near infrared radiation (NIR) and then summed to yield total intercepted radiation.Down-scattered radiation (Q sc ) is the difference between total and direct beam: Finally, the mean flux densities on sunlit (Q sl ), and shaded (Q sh ), leaves, respectively, are calculated as: [36] TREES was modified to explicitly account for spatial variabity in incoming solar radiation at each 30 min time step to incorporate the effects of canopy heterogeneity on incident radiation.Transmittance for diffuse radiation is integrated over all zenith angles and because of this canopy heterogeneity contributes less to variability in diffuse radiation, as such the impacts of this heterogeneity on diffuse radiation were considered negligible and were not accounted for.Incident beam radiation, Q ob , as described above, was multiplied by the proportion of light transmitted through the canopy t bc , which can be defined as follows [Campbell and Norman, 1998]: where m is leaf area density (m 2 m −3 ) and S(Y;) is the path length of light rays through the canopy.S(Y;) was estimated by plotting the Sun's path on hemispherical outputs of path length for each tree generated using tRAYci (Figure 3).
[37] Acknowledgments.Funding for this study was from the National Science Foundation, Hydrologic Sciences Program (EAR-0405306 to D.S.M., EAR-0405381 to B.E.E., and EAR-0405318 to E.L.K.).D.S.M. also acknowledges funding from the Department of Energy (DOE) Office of Biological and Environmental Research, National Institute for Climate Change Research (NICCR) Midwestern region subagreement 050516Z20.
In addition, M.M.L. was supported by the NSF IGERT program.The statements made in this manuscript reflect the views of the authors and do not necessary reflect the views of NSF.We are also grateful to personnel at Kemp Natural Resources Station, and to Dave Roberts for their assistance with fieldwork.We would like to thank two anonymous reviewers for helpful comments on earlier versions of the manuscript.

Figure 1 .
Figure 1.A conceptual diagram illustrating the models compared in this study.In the big leaf model (a) TREES requires a combination of micrometeorological and environmental inputs, and ecophysiological parameters to estimate fluxes of carbon and water.Unknown parameters G Sref and m are repeatedly sampled from a predefined parameter space (n = 10,000).E C from each of the resulting 10,000 simulations is compared to observed E C in order to identify the best simulation to use for analysis.In the light model (b) inputs from stand inventories are used to calculate the distance a ray of sunlight must travel through the neighboring canopy to reach the crown of interest using the tRAYci model.The TREES model is modified to accept path length as an input, which is used to constrain incident radiation and partition LAI into sunlit and shaded portions.All other inputs and outputs for TREES are the same, and an identical model evaluation procedure are followed to identify the best out of 10,000 simulations to be used for analysis.Note that unique estimates of LAI for each tree, based on allometry, are used as model inputs.

Figure 2 .
Figure 2. Diurnal plots of transpiration from three representative individuals covering the range of observed rates, and of D and Q 0 (a and b) for Aspen, and (c and d) Maple sites.

Figure 3 .
Figure 3. Simulated E C plotted against observed E C for all aspen and maple individuals respectively.Dashed lines represent 1:1 relationship and solid lines are linear regressions.Black lines indicate big leaf model regressions and gray lines indicate light model regressions.big leaf model r 2 = 0.72, 0.84 and slopes of 0.90 and 0.94 for aspen and maple respectively.Light model r 2 = 0.75, 0.80 and slopes of 0.90 and 0.91 for aspen and maple respectively.

aFigure 4 .
Figure 4. Diurnal plots of observed and simulated E C and simulated photosynthetic light response curves for maple trees (a and b) J11, (c and d) E6, and (e and f) A5.

Figure 5 .
Figure 5.Diurnal plots of observed and simulated E C and simulated photosynthetic light response curves for aspen trees H6 (a) and (b), A11 (c) and (d), and D7 (e) and (f).

Table 1 .
Summary of Stand Inventory Data Used to Model Light and Assess Competition a Species accounting for less than 1% of the sampled population are not listed.bMeasuredat Aspen site, estimated at maple using available allometry.cTotal count includes species not listed.d Trees with multiple stems that split below 1.3 m are counted separately.

Table 2 .
Mean Parameter Values