Dynamic Simulation of the Crown Net Photosynthetic Rate for Young Larix olgensis Henry Trees

Numerical integration of the instantaneous net photosynthetic rate (An) is a common method for calculating the long-term CO2 uptake of trees, and accurate dynamic simulation of the crown An has been receiving substantial attention. Tree characteristics are challenging to assess given their aerodynamically coarse crown properties, spatiotemporal variation in leaf functional traits and microenvironments. Therefore, the variables associated with the dynamic variations in the crown An must be identified. The relationships of leaf temperature (Tleaf), the vapor pressure deficit (VPD), leaf mass per area (LMA) and the relative depth into the crown (RDINC) with the parameters of the photosynthetic light-response (PLR) model of Larix olgensis Henry were analyzed. The LMA, RDINC and VPD were highly correlated with the maximum net photosynthetic rate (Amax). The VPD was the key variable that mainly determined the variation in the apparent quantum yield (AQY). Tleaf exhibited a significant exponential correlation with the dark respiration rate (Rd). According to the above correlations, the crown PLR model of L. olgensis trees was constructed by linking VPD, LMA and RDINC to the original PLR equation. The model performed well, with a high coefficient of determination (R2) value (0.883) and low root mean square error (RMSE) value (1.440 μmol m−2 s−1). The extinction coefficient (k) of different pseudowhorls within a crown was calculated by the Beer–Lambert equation based on the observed photosynthetically active radiation (PAR) distribution. The results showed that k was not a constant value but varied with the RDINC, solar elevation angle (ψ) and cumulative leaf area of the whole crown (CLA). Thus, we constructed a k model by reparameterizing the power function of RDINC with the ψ and CLA, and the PAR distribution within a crown was therefore well estimated (R2 = 0.698 and RMSE = 174.4 μmol m−2 s−1). Dynamic simulation of the crown An for L. olgensis trees was achieved by combining the crown PLR model and dynamic PAR distribution model. Although the models showed some weakened physiological biochemical processes during photosynthesis, they enabled the estimation of long-term CO2 uptake for an L. olgensis plantation, and the results could be easily fitted to gas-exchange measurements.


Introduction
Photosynthesis is one of the most important elements in many ecophysiological models [1][2][3] as photosynthesis plays a vital role in the material cycle and energy flow of forest ecosystems [1,4,5].Photosynthesis supports metabolism in leaves and provides nutrition to stems and other plant parts for growth [6].For an individual tree, a good understanding of crown photosynthesis also informs physiological principles for forest management, such as artificial pruning for young forests [7], which can improve the economic value of stems [8,9].However, the simulation of photosynthesis is generally challenging due to the spatiotemporal variability in leaf functional traits [6,7,[10][11][12][13][14][15][16][17], microenvironments [5,[18][19][20], foliage ages [12,[20][21][22]] and branch factors [23].Many mechanistic models, including the big-leaf model [24], the multistory model [25] and the two-leaf model [26], which were developed from Farquhar's scheme [27], accurately describe the biochemical process of crown photosynthesis.These models have been widely applied in many ecophysiological process-based models [1,28].However, these models are usually complex as they require empirical calculations of species-specific photo-and biochemical parameters that are subject to environmental control [29], and the behaviors of the parameters are not easily predicted from their formulations [3,20].Semiempirical models, such as the photosynthetic light-response (PLR) curve model, serve as an intermediate option to define photosynthesis according to a mathematical description of a single biochemical process between the net photosynthetic rate (A n ) and photosynthetically active radiation (PAR) [29].
PLR curves vary considerably among species and among individuals of the same species growing in different environments [7,30].To expand the application of the PLR model from a single leaf to larger scales, many researchers have tried to link leaf functional traits and environmental conditions to the main photosynthetic indicators [11,16,20,29,31,32].Leaf mass per area (LMA) and its reciprocal variable, the specific leaf area (SLA), are both important leaf functional traits that have been suggested to be closely associated with mesophyll conductivity to CO 2 [33,34], consequently impacting photosynthetic indicators [13,[35][36][37].Thus, some researchers have estimated photosynthetic indicators using LMA according to their allometric relationships [11,16,38,39], but photosynthetic indicators can be more accurately predicted if the light condition is included [40].The nitrogen (N) content is one of the most important factors related to leaf photosynthetic capability [11,12,16,41].Xu et al. [32] used a chlorophyll meter (SPAD) to incorporate a popular indicator reflecting N into the PLR model to overcome the large coefficient variations between individual leaves with different N statuses.Photosynthesis is also sensitive to changes in surrounding environmental conditions [42][43][44][45].Ignoring environmental effects may reduce the precision and applicability of large-scale models [16].Calama et al. [20] and Mayoral et al. [29] successfully expanded the application scope of the PLR model by linking leaf temperature (T leaf ) and the global site factor (GSF) to the parameters.Photosynthetic indicators within crowns also vary with the vertical position [38,[46][47][48][49][50]; therefore, simulation of crown photosynthesis should comprehensively consider leaf functional traits, environmental conditions and vertical position within a crown.Liu et al. [7] estimated the response of the A n to PAR at different crown positions and under different environmental conditions by modifying the PLR model via reparameterization according to the relationships between the main parameters and LMA, T leaf , the vapor pressure deficit (VPD), and the relative depth into the crown (RDINC).
Light transmission also plays a crucial role in calculating CO 2 uptake due to the sensitive response of the photosynthetic process to light.Thus, light transmission has been simulated in many forest dynamics models, e.g., FORSKA [51]; FORMIND [52]; LPJ [53]; PHOTOS [1]; FAREST [54]; JULES [55].Beer's law was the core theory in these models, where the Beer-Lambert equation was the basic equation [56][57][58].The extinction coefficient (k) and leaf area index (LAI) are the only two essential factors that depend on the distributions of the leaf angle and leaf area, respectively.Campbell [59] used the ellipsoidal leaf angle density function to yield good approximations of leaf angle distributions in real plant canopies and calculated the formula of k.Then, Campbell [60] and Wang and Jarvis [61] developed empirical equations for k involving the mean leaf angle.The LAI was also modified by using the clumping index [62] because leaves are generally non-randomly distributed in a canopy or even in an individual crown [63].Because leaf angle distributions usually exhibit high spatial and temporal variability [64], k is not unique [1,63,65], although some previous studies directly use a constant k value [14,[52][53][54][55]66].In addition, light is also intercepted by other phytoelements, such as twigs and branches [61,62], which hinders an accurate simulation of light transmission.
Larix olgensis is one of the main tree species used for afforestation in northeastern China, and this species accounts for 36% and 37% of the total area and total volume, respectively, of all plantations located in this region.A good understanding of the dynamic A n of L. olgensis trees is significant for studying carbon storage and the carbon cycle in the main forest type in this area.Annual or long-term CO 2 uptake is generally determined by integrating the instantaneous A n [1,14,28]; thus, an accurate dynamic simulation of the A n is crucial.The objectives of this investigation are (1) to expand the PLR model from a single leaf to a tree crown with varied leaf functional traits and environmental conditions for young planted L. olgensis trees, (2) to model the radiation transfer within a crown, and (3) to simulate the dynamic A n of a multilayer crown through the growing season.

Site Description and Sample Tree Selection
The experiments were conducted at the Maoershan Forest Farm (127 • 18 0" E-127 • 41 6" E and 45 • 2 20" N-45 • 18 16" N) in northeastern China.The mean annual air temperature, the mean annual precipitation, and the mean elevation in this region are 2.0 • C, 525 mm, and 400 m, respectively.The soil type is Eutroboralf, and the total forest coverage is approximately 83.3% of the area, 14.7% of which consists of plantations.
Five sample trees, each with a similar diameter at breast height (DBH) to quadratic mean diameter (Dg), representing the average state of each sample plot (20 m × 30 m) were selected according to the measurements of basic tree factors for each tree.Scaffoldings were constructed around each sample tree, needle samples for gas-exchange measurements were selected at different positions within the crown, and this process is described in detail in Liu and Li [17].

Measurements of the PLR Curves
The PLR curves were generated using fully expanded healthy needles and a portable photosynthesis system (LI-6400XT, LI-COR, Inc., Lincoln, NE, USA) coupled with a standard light-emitting diode (LED) light source (Li-6400-02B) at 17 PAR levels (2200, 2000, 1800, 1600; 1400, 1200, 1000, 800, 600, 400, 200, 150, 100, 60, 40, 15 and 0 µmol m −2 s −1 ).Sample needle clusters were acclimated for 20 min at a CO 2 concentration of 370 µmol m −2 s −1 , which approached the actual concentration around the needles, although this concentration was slightly lower than the ambient conditions, by using a CO 2 mixer (Li-6400-01) to maintain a stable CO 2 supply in the chamber.In our previous study, we proved that A n in different vertical positions within crowns of Larix olgensis trees tended to be stable at a PAR of 1400 µmol m −2 s −1 [67]; thus, 1400 µmol m −2 s −1 was used to induce all needle clusters for at least 10 min to ensure that the potential photosynthetic capacities of all needle clusters could be activated.The needle clusters were allowed to equilibrate for a minimum of 2 min at each step before the data were logged.The data for the PLR curves were measured every half month during the growing season (approximately from 15th May to 10th September) in 2016 and 2017.The measurements were conducted from 7:00 to 17:00.

Measurements of the RDINC and LMA
The depth into the crown (DINC) of the needle clusters and the data for the PLR curves were measured.Then, the RDINC of each needle cluster was calculated by dividing the DINC by crown depth (the distance from the tree tip to the base of its live crown, CD).Once the photosynthetic gas-exchange measurements were completed, corresponding needle clusters were scanned immediately and then surveyed with image analysis software (Image-Pro Plus 6.0, Media Cybernetics, Bethesda, MD, USA) in the laboratory, resulting in a projected leaf area.Then, immediately after their collection in the field, the corresponding needles were dried to a constant weight at 85 • C and weighed.LMA was calculated by dividing the needle dry mass by the corresponding projected leaf area.The statistics for the measurement data for the PLR curves and the corresponding T leaf , VPD and LMA are listed in Table 1.

Measurements of PAR within the Crown
PAR measurements were conducted after the leaves fully expanded.First, the incident PAR at the top of the crown was measured using a PAR sensor on the handheld part of the photosynthesis system (LI-6400XT, LI-COR, Inc., Lincoln, NE, USA).Then, the PAR in each pseudowhorl was measured in four directions depending on the branch azimuth (ϕ) (0 ≤ ϕ < 45 and 315 ≤ ϕ < 360 defined as north; 45 ≤ ϕ < 135 defined as east; 135 ≤ ϕ < 225 defined as south; and 225 ≤ ϕ < 315 defined as west) with 10 records every 30 cm in horizontal distance from the tree stem to the outermost part of the pseudowhorl.Each measurement for the whole tree crown was conducted in at 8:00, 10:00, 12:00, 14:00 and 16:00 on a clear and cloudless day for each sample tree.The solar elevation angle (ψ) at each measurement was also recorded.The data in this section were used to fit the k model.

Calculation of the Leaf Area and LAI
At the end of the growing season (early September), three sample trees were cut down for measurement.The CD of each sampled tree was divided into several segments based on the pseudowhorls [68] from top to bottom.The branch diameter (BD), branch length (BH), branch chord length (BC), ϕ, branch angle (θ) and DINC of all branches in each pseudowhorl were measured.Then, all branches were cut and weighed.Next, in each pseudowhorl, one branch with the average weight was selected as the representative branch, and the branch and needle clusters were separated and weighed.The ratio of needle-cluster weight to branch weight (RNB) in each segment was calculated.The needle clusters were then sampled (approximately 50-100 g) randomly, weighed and taken to the laboratory for dry matter proportion (DMP) determination.The LMA for each corresponding pseudowhorl was also obtained by additional needle-cluster samples (approximately 10 groups and each group was 0.1 g).Finally, the LA in each pseudowhorl could be calculated as follows: where TW is the total weight of the branches and needle clusters in each pseudowhorl, and i is the serial number of the pseudowhorl.The relative crown projected area (RCA) of each pseudowhorl was calculated as a circle using the mean projected length of the four longest branches in four directions.The relative cumulative LAI (RCLAI) of each pseudowhorl was calculated by dividing the relative cumulative leaf area (RCLA) by the RCA.The statistics of the measured PAR, LA and RCLA are listed in Table 2.
Table 2. Summary of the main measured variables.The following variables are given: the photosynthetically active radiation in each pseudowhorl (PAR i ), incident PAR at the top of the crown (PAR t ), the solar elevation angle (ψ), leaf area in each pseudowhorl (LA i ), relative cumulative leaf area index in each pseudowhorl (RCLAI i ).i is the serial number of pseudowhorls belonging to five planted Larix olgensis trees.

Meteorological Data Collection
The instantaneous T leaf , air temperature (T air ) and VPD were recorded simultaneously when the PLR curves were measured using the photosynthesis system (LI-6400XT, LI-COR, Inc., Lincoln, NE, USA).The dynamic data of T air , the VPD and PAR (recorded every 30 min each day) were obtained from the National Science and Technology Infrastructure Program (http://mef.cern.ac.cn/).The data in this section were used to simulate the dynamic crown A n of L. olgensis.

Model Descriptions
In this study, the modified Mitscherlich equation [69] was used to model the PLR curve: where A max is the maximum net photosynthetic rate (µmol m −2 s −1 ), AQY is the apparent quantum yield, and R d is the dark respiration rate (µmol m −2 s −1 ).
According to a previous study [7], the correlations of the parameters (the A max , AQY and R d ) with the environmental conditions (T leaf and the VPD), leaf traits (LMA) and spatial position (the RDINC) were analyzed, and the original PLR curve (Equation ( 2)) was modified to the following form: where f (x), g(x) and h(x) represent the functions of T leaf , VPD, LMA, and RDINC or their combinations.
The Beer-Lambert equation was used to model the PAR distribution within the crown: where i is the serial number of the pseudowhorl, PARt is the incident PAR at the top of the crown, and ψ is the solar elevation angle.The k value in each pseudowhorl was calculated with Equation (4) based on the observed PAR i , PAR t , RCLAI i and ψ.Then, the vertical patterns of the k value for the sample trees were modelled based on the relationship between k and the RDINC, ψ and CLA as shown below: where f (RDINC) represents the function of the RDINC.

Model Assessment and Validation
The crown PLR model was fitted to the entire data set, and the models were validated using the jackknife technique, in which a PLR model was constructed using four trees and the other one tree was used to validate the model.We repeated the process for five times to make sure that all the five sample trees have been used as validation data at least once.The model assessment was based on the coefficient of determination (R 2 ) and the root mean square error (RMSE) (Equations ( 6) and ( 7)), and model validation was evaluated based on R 2 , RMSE, the mean error (ME) and precision estimation (P) (Equations ( 8) and ( 9)) as follows: where y i is the observed value; y is the mean of the observed value; ŷi is the predicted value; n is the number of observations; and p is the number of parameters.

Correlations between Photosynthetic Parameters and the RDINC, T leaf , VPD and LMA
Figure 1 shows the relationships between the fitted parameters (the A max , AQY and R d ) for the Mitscherlich equation and the RDINC, T leaf , VPD and LMA in the 752 needle clusters from the five sample trees.The A max exhibited strong correlations with the RDINC and LMA, in which A max was negatively correlated with the RDINC but positively correlated with LMA.Although the relationships between the A max and environmental conditions (T leaf and the VPD) were also significant, their correlations were weak.The AQY exhibited a strong negative correlation with the VPD but weak correlations with T leaf .The R d was positively correlated with T leaf but negatively correlated with the RDINC.In addition, the VPD and LMA were also positively correlated with the R d .
where yi is the observed value;  is the mean of the observed value;  is the predicted value; n is the number of observations; and p is the number of parameters.

Correlations between Photosynthetic Parameters and the RDINC, Tleaf, VPD and LMA
Figure 1 shows the relationships between the fitted parameters (the Amax, AQY and Rd) for the Mitscherlich equation and the RDINC, Tleaf, VPD and LMA in the 752 needle clusters from the five sample trees.The Amax exhibited strong correlations with the RDINC and LMA, in which Amax was negatively correlated with the RDINC but positively correlated with LMA.Although the relationships between the Amax and environmental conditions (Tleaf and the VPD) were also significant, their correlations were weak.The AQY exhibited a strong negative correlation with the VPD but weak correlations with Tleaf.The Rd was positively correlated with Tleaf but negatively correlated with the RDINC.In addition, the VPD and LMA were also positively correlated with the Rd.

Simulation of the Responses of the A n to PAR under Different Conditions
The PLR curve model (Equation ( 3)) was determined according to the relationships between the fitted parameters (the A max , AQY and R d ) for the PLR curve equation and the RDINC, T leaf , VPD and LMA (Figure 1).To make sure that the scaled parameters are consistent with the ones retrieved at leaf level by fitting the PLR curves, the final crown PLR curve model was determined as follow: where a 0 , a 1 , b 0 , b 1 and R d are unknown parameters to be estimated.Table 3 shows the parameter estimates for the nonlinear least-squares fit of model 10 and the goodness of fit and validation statistics for the model.All the parameters were significant at a 0.01 level, and the model presented satisfactory model assessment in modeling the crown PLR curves, with a high

Vertical Patterns of PAR and k within the Crown
Figure 2 shows the vertical patterns of light transmittance (the ratio of PAR in the crown to the incident PAR, LT) for five planted L. olgensis trees.Although the LT decreased with increases in the RDINC in the five sample trees, the patterns still showed minimal differences among different trees, which was probably associated with the different PAR t and RCLAI in each pseudowhorl (Table 2).This scenario also caused a vertical variation in k, where k increased with an increasing RDINC (Figure 3).In the top of the crown, the k for the five sample trees was relatively stable at approximately 0.2.Compared to the k values in the top of the crown, the k values in the bottom of the crown k values in tree I, tree II, tree III, tree IV and tree V were 1.9-, 1.8-, 4.8-, 16.1-and 7.2-times higher, respectively.In addition, k was also different among the five sample trees, as the mean k values were 0.39, 0.32, 0.51, 0.61 and 0.58 in tree I, tree II, tree III, tree IV and tree V, respectively.

Estimation of k
According to the nonlinear correlation between k and the RDINC, we modeled k for the five sample trees based on the power function as  =  • RDINC .Then, the parameters (a0 and b0) were reparameterized by linking ψ and CLA to model k.The k model was finally constructed as follows: where d0, d1, d2, e0 and e1 are unknown parameters to be estimated.
Table 4 shows the parameter estimates for the nonlinear least-squares fit of model 11 and the goodness of fit statistics for the model.All the parameters were significant at a 0.01 level, and the model performed well in modeling the multiwhorl k, with a high R 2 value (0.569) and a low RMSE

Estimation of k
According to the nonlinear correlation between k and the RDINC, we modeled k for the five sample trees based on the power function as k = a 0 •RDINC b 0 .Then, the parameters (a 0 and b 0 ) were reparameterized by linking ψ and CLA to model k.The k model was finally constructed as follows: where d 0 , d 1 , d 2 , e 0 and e 1 are unknown parameters to be estimated.Table 4 shows the parameter estimates for the nonlinear least-squares fit of model 11 and the goodness of fit statistics for the model.All the parameters were significant at a 0.01 level, and the model performed well in modeling the multiwhorl k, with a high R 2 value (0.569) and a low RMSE value (0.159).Meanwhile it also performed well in validation, with high values of R 2 (0.569) and P% (93.15) and low values of RMSE (0.170) and ME (0.038).

Dynamic Simulation of the A n within the Crown throughout the Whole Growing Season
With the dynamic meteorological data, the PAR in each pseudowhorl within the crown was calculated based on the multiwhorl k model, and consequently, the dynamic A n per area (A n-area ) in Forests 2019, 10, 321 each pseudowhorl was simulated according to model 10. Figure 4a shows the dynamic A n-area at different pseudowhorls during July (tree I was the example), which clearly describe the temporal variation in the A n-area in multiwhorls within a crown at a diurnal scale and monthly scale and even at a seasonal scale.The A n-area decreased with an increasing RDINC, but the A n per pseudowhorl (A n-pw ) was greatest in the middle crown (Figure 4b).We calculated the daily net photosynthetic production (NPP D ) of the whole crown according to the dynamic simulation of the crown A n based on Equation ( 8) (NPP D -I) and Equation ( 2) (NPP D -II) respectively (Figure 5).The NPP D that was calculated based on our modified PLR model (NPP D -I) was highly consistent with the dynamic changes in the cumulative daily PAR (PAR D ), which conformed to the natural reality.However, the NPP D that was calculated based on the original PLR model (NPP D -II) was weakly correlated with the PAR D and was obviously an overestimation.The NPP D exhibited obvious fluctuations during the growing season but was highly consistent with the dynamic changes in the cumulative daily PAR (PAR D ).

Dynamic Simulation of the An within the Crown throughout the Whole Growing Season
With the dynamic meteorological data, the PAR in each pseudowhorl within the crown was calculated based on the multiwhorl k model, and consequently, the dynamic An per area (An-area) in each pseudowhorl was simulated according to model 10. Figure 4(a) shows the dynamic An-area at different pseudowhorls during July (tree I was the example), which clearly describe the temporal variation in the An-area in multiwhorls within a crown at a diurnal scale and monthly scale and even at a seasonal scale.The An-area decreased with an increasing RDINC, but the An per pseudowhorl (Anpw) was greatest in the middle crown (Figure 4b).We calculated the daily net photosynthetic production (NPPD) of the whole crown according to the dynamic simulation of the crown An based on Equation ( 8) (NPPD-I) and Equation (2) (NPPD-II) respectively (Figure 5).The NPPD that was calculated based on our modified PLR model (NPPD-I) was highly consistent with the dynamic changes in the cumulative daily PAR (PARD), which conformed to the natural reality.However, the NPPD that was calculated based on the original PLR model (NPPD-II) was weakly correlated with the PARD and was obviously an overestimation.The NPPD exhibited obvious fluctuations during the growing season but was highly consistent with the dynamic changes in the cumulative daily PAR (PARD).

Discussion
For trees, the photosynthetic properties of foliage exhibit vertical variation within a crown [47,50,70], and this phenomenon is intimately associated with leaf traits [10,15,23] and the surrounding microenvironment [18].Our results showed that the Amax and Rd were negatively

Discussion
For trees, the photosynthetic properties of foliage exhibit vertical variation within a crown [47,50,70], and this phenomenon is intimately associated with leaf traits [10,15,23] and the surrounding microenvironment [18].Our results showed that the A max and R d were negatively correlated with the RDINC (Figure 1), which is consistent with other studies [35,48,49,[71][72][73][74].This phenomenon is consistent with the theory of low light acclimation [23,70].LMA generally demonstrated a strong correlation with photosynthesis [31,[35][36][37] due to its close association with mesophyll conductivity to CO 2 [33,34], and LMA has also served as an indicator of the potential growth rate [75].Figure 1 shows that LMA exhibited a positive correlation with the A max , which is consistent with previous studies [38,[76][77][78].In contrast, LMA is usually negatively correlated with the A max in large trees [34,[79][80][81].However, the LMA in large trees can be extremely high, which may hinder mesophyll conductivity to CO 2 and consequently reduce the A max [34,35].T leaf has been extensively reported to be parabolically correlated with the A max [20,29,[82][83][84][85], and the optimum temperature (Topt) for the A max was in the range of 15-30 • C for most temperate species [12].Our results (Figure 1) are also in agreement with the results from the above reports; however, the correlation between T leaf and the A max was weak, and the optimum temperature for A max was approximately 30 • C, which is consistent with the optimum temperature for the A max of Pinus pinea [20].The R d demonstrated a significant exponential correlation with T leaf (Figure 1), which is also in line with the results of other previous studies [5,49,84,86].Our results showed that the AQY was strongly and negatively correlated with the VPD (Figure 1), presumably due to the limitation that a high VPD imposes on stomatal conductance [19,[87][88][89].Similar to Cannell and Thornley [90], we also found a decreasing relationship between the AQY and T leaf , although some previous studies have suggested that the AQY is independent of temperature [20].
The A max , AQY and R d are not only important photosynthetic indicators but also core parameters for determining the shape of the PLR curve (Equation ( 2)).Thus, many studies have linked leaf traits [11,16,32,91] or environmental conditions [20,29] in the PLR models based on their correlations with the A max , AQY and R d to resolve the problem of parameter variation among different species or environments.In this study, we comprehensively considered the influence of spatial position, leaf traits and environmental conditions on leaf photosynthesis and then linked the RDINC, LMA and VPD in the PLR model to extend the PLR model from a single leaf to the whole crown and render the model sufficiently flexible to be applied in different conditions.The Beer-Lambert equation was the basic equation [56][57][58] used to model PAR transmission within the crown, in which the k value played a vital role.However, the k value was challenging to determine in the trees due to their aerodynamically coarse crown properties.The ellipsoidal leaf angle density function is generally used to estimate the leaf angle distributions in real plant crowns to calculate k values [59].Campbell [60] and Wang and Jarvis [61] developed empirical equations for k involving the mean leaf angle, but the equations are time-consuming and complex in terms of leaf angle measurement.Thus, many studies directly use a constant k value of 0.4 [54], 0.5 [14,55], 0.52 [66] or 0.7 [52] to simulate PAR transmission within crowns.Our results showed that the k value was not constant (Figure 3), which was also observed in other studies [63,65].In addition, we modeled k based on its exponential correlation with the RDINC for each sample tree and further constructed the k model (Equation ( 11)) by reparameterization using ψ and CLA. Figure 6 shows the performances of different constant k values in terms of predicting PAR within the crowns of the five sample trees.The k model performed better when k was assigned a value of 0.33, with a greater R 2 (0.37) and smaller RMSE (252.31µmol m −2 s −1 ).We contrasted the results of the estimated PAR using different methods for determining k (Figure 7), and the results showed that our k model performed the best among these models (Figure 7a).In addition, the k value was unsatisfactorily determined based on the leaf angle (Figure 7c,d) for L. olgensis trees.
values in terms of predicting PAR within the crowns of the five sample trees.The k model performed better when k was assigned a value of 0.33, with a greater R 2 (0.37) and smaller RMSE (252.31μmol m −2 s −1 ).We contrasted the results of the estimated PAR using different methods for determining k (Figure 7), and the results showed that our k model performed the best among these models (Figure 7a).In addition, the k value was unsatisfactorily determined based on the leaf angle (Figure 7c,d) for L. olgensis trees.The crown An is difficult to simulate, specifically in trees with their aerodynamically coarse microenvironments and seasonal and spatial changes in leaf traits.Thus, a crown An model should be sufficiently enough to consider account the diurnal, seasonal and spatial variations in the An.Previous studies have suggested that incorrect simulations usually occur when modeling long-term carbon uptake if the spatial and seasonal variations in photosynthesis [46,[92][93][94][95] or leaf traits [96,97] are neglected.Some scientists have successfully incorporated leaf traits [11,16,29,32,91] and environmental conditions [20,29] into the PLR model to extend the model from a single leaf to a larger scale.However, these models are limited in their application because they only describe the response of the An to PAR in different conditions and do not describe the dynamic simulation of the An.In this study, we constructed a k model to simulate the dynamic PAR within a crown to achieve a dynamic simulation of the An based on dynamic meteorological data (Figure 5).
For a clear understanding of the influence of the k value determined based on different methods on the NPPD, we compared the NPPD results that were calculated based on two methods (Figure 8) for tree I.The total NPPD that was calculated by method II was overestimated by 14.1% compared to The crown A n is difficult to simulate, specifically in trees with their aerodynamically coarse microenvironments and seasonal and spatial changes in leaf traits.Thus, a crown A n model should be sufficiently enough to consider account the diurnal, seasonal and spatial variations in the A n .Previous studies have suggested that incorrect simulations usually occur when modeling long-term carbon uptake if the spatial and seasonal variations in photosynthesis [46,[92][93][94][95] or leaf traits [96,97] are neglected.Some scientists have successfully incorporated leaf traits [11,16,29,32,91] and environmental conditions [20,29] into the PLR model to extend the model from a single leaf to a larger scale.However, these models are limited in their application because they only describe the response of the A n to PAR in different conditions and do not describe the dynamic simulation of the A n .In this study, we constructed a k model to simulate the dynamic PAR within a crown to achieve a dynamic simulation of the A n based on dynamic meteorological data (Figure 5).
For a clear understanding of the influence of the k value determined based on different methods on the NPP D , we compared the NPP D results that were calculated based on two methods (Figure 8) for tree I.The total NPP D that was calculated by method II was overestimated by 14.1% compared to the total NPP D that was calculated by method I.The NPP D results for each pseudowhorl showed that method I underestimated the value in the upper pseudowhorl compared to method II but overestimated the value in the lower pseudowhorl compared to method II, indicating that the NPP D calculated by method II may have yielded an incorrect result regarding the contribution of photosynthate from branches in pseudowhorls to trunks.Consequently, this method resulted in an incorrect judgment for the effective crown (i.e., the basis for the scientific, artificial pruning of young forests).In summary, the results of this study show that environmental conditions (Tleaf and the VPD), leaf traits (LMA) and spatial positions (the RDINC) are the key variables causing the dynamic variation in crown PLR curves for L. olgensis trees.Thus, Tleaf, the VPD, LMA and the RDINC should to be linked to the parameters of the PLR model to construct a crown PLR model of L. olgensis trees.The k value exhibits obvious vertical variation within a crown, indicating that a constant value of k will result in incorrect PAR estimations.Therefore, a k model was established to simulate the dynamic PAR within a crown based on the Beer-Lambert function.Finally, the dynamic simulation of the crown An for L. olgensis trees was conducted by combining the crown PLR model and dynamic PAR distribution model.Compared to the mechanistic photosynthesis models [24-27] that have been widely used in forest dynamics models [51][52][53][54][55], our model weakened the descriptions of some physiological and biochemical processes, rendering the model more convenient for application than the mechanistic models.The model can be further improved based on data from trees of different sizes and can serve as a foundation for estimating the long-term CO2 uptake in L. olgensis plantations.

Conclusions
The PLR curves and PAR were measured within multilayer tree crowns of five planted L. olgensis trees during the whole growing season.The crown PLR model of L. In summary, the results of this study show that environmental conditions (T leaf and the VPD), leaf traits (LMA) and spatial positions (the RDINC) are the key variables causing the dynamic variation in crown PLR curves for L. olgensis trees.Thus, T leaf , the VPD, LMA and the RDINC should to be linked to the parameters of the PLR model to construct a crown PLR model of L. olgensis trees.The k value exhibits obvious vertical variation within a crown, indicating that a constant value of k will result in incorrect PAR estimations.Therefore, a k model was established to simulate the dynamic PAR within a crown based on the Beer-Lambert function.Finally, the dynamic simulation of the crown A n for L. olgensis trees was conducted by combining the crown PLR model and dynamic PAR distribution model.Compared to the mechanistic photosynthesis models [24-27] that have been widely used in forest dynamics models [51][52][53][54][55], our model weakened the descriptions of some physiological and biochemical processes, rendering the model more convenient for application than the mechanistic models.The model can be further improved based on data from trees of different sizes and can serve as a foundation for estimating the long-term CO 2 uptake in L. olgensis plantations.

Conclusions
The PLR curves and PAR were measured within multilayer tree crowns of five planted L. olgensis trees during the whole growing season.The crown PLR model of L. olgensis trees was constructed by linking VPD, LMA and RDINC to the original PLR equation.The k model was constructed based on the ψ, CLA and RDINC, and consequently the dynamic PAR distribution can be simulated.Finally, the dynamic simulation of multilayer crown A n for L. olgensis trees was conducted by combining the crown PLR model and dynamic PAR simulation.Compared with the traditional mechanism models, our model weakened other physiological biochemical processes but with simpler model structure, fewer variables to be measured and guaranteed well estimation accuracy.Based on our model the NPP in each pseudowhorl through the growing season can be calculated by numerical integration of the dynamic A n , thus serving as a foundation for determining the ratio of the functional crown, which could provide a scientific and effective guidance in artificial pruning for young forests.

Figure 1 .
Figure 1.Relationships between the fitted parameters (The maximum net photosynthetic rate, A max , µmol m −2 s −1 ; the apparent quantum yield, AQY; and the dark respiration rate, R d , µmol m −2 s −1 ) for the modified Mitscherlich equation and the relative depth into the crown (RDINC), leaf temperature (T leaf , • C), vaper pressure deficit (VPD, kPa) and leaf mass per area (LMA, g m −2 ) in the 752 needle clusters from the five sample trees (open circles).Color circles containing values show Spearman's correlation coefficients (α = 0.05).

Forests 2018, 9 , 20 Figure 2 .
Figure 2. Vertical patterns of light transmittance (The ratio of photosynthetically active radiation in the crown to the incident photosynthetically active radiation, LT) for the planted Larix olgensis Henry trees.(a~e) represents Tree I ~ Tree V.The dark circles are the mean LT in each pseudowhorl, and the black bars are standard deviations.

Figure 2 .Forests 2019, 10 , 321 9 of 18 Figure 2 .
Figure 2. Vertical patterns of light transmittance (The ratio of photosynthetically active radiation in the crown to the incident photosynthetically active radiation, LT) for the planted Larix olgensis Henry trees.(a~e) represents Tree I ~Tree V.The dark circles are the mean LT in each pseudowhorl, and the black bars are standard deviations.

Figure 3 .
Figure 3. Vertical patterns of the extinction coefficient (k) for the five planted Larix olgensis trees.(a~e) represents Tree I ~ Tree V.The RDINC is the relative depth into the crown; the dark circles are mean k in each pseudowhorl; the black bars are standard deviations; and the dotted lines are the trendlines following the power function.

Figure 3 .
Figure 3. Vertical patterns of the extinction coefficient (k) for the five planted Larix olgensis trees.(a~e) represents Tree I ~Tree V.The RDINC is the relative depth into the crown; the dark circles are mean k in each pseudowhorl; the black bars are standard deviations; and the dotted lines are the trendlines following the power function.

Figure 4 .
Figure 4.The dynamic simulation of (a) the net photosynthetic rate (A n ) per area (A n-area ) and (b) the total A n per pseudowhorl (A n-pw ) (only tree I was used as an example during July).Three pseudowhorls with relative depths into the crown (RDINC) of 0.22, 0.49 and 0.88, representing the upper, middle and lower crowns, respectively, were selected.

Figure 4 .
Figure 4.The dynamic simulation of (a) the net photosynthetic rate (An) per area (An-area) and (b) the total An per pseudowhorl (An-pw) (only tree I was used as an example during July).Three pseudowhorls with relative depths into the crown (RDINC) of 0.22, 0.49 and 0.88, representing the upper, middle and lower crowns, respectively, were selected.

Figure 5 .
Figure 5.The daily net photosynthetic production of the whole crown (NPPD) based on Equation (8) (NPPD-I) and Equation (2) (NPPD-II) for tree I and the corresponding cumulative daily photosynthetically active radiation (PARD).

Figure 5 .
Figure 5.The daily net photosynthetic production of the whole crown (NPP D ) based on Equation (8) (NPP D -I) and Equation (2) (NPP D -II) for tree I and the corresponding cumulative daily photosynthetically active radiation (PAR D ).

Figure 6 .
Figure 6.The performance of the estimation of the photosynthetically active radiation (PAR) distribution within crowns with different extinction coefficient (k) values.R 2 is the coefficient of determination, and RMSE is the root mean square error.

Figure 6 .
Figure 6.The performance of the estimation of the photosynthetically active radiation (PAR) distribution within crowns with different extinction coefficient (k) values.R 2 is the coefficient of determination, and RMSE is the root mean square error.

Forests 20 Figure 7 .
Figure 7.Comparison of the observed and simulated photosynthetically active radiation (PAR) within a crown: (a) PAR is simulated with the k model that was constructed in this study; (b) PAR is simulated with k values ranging from 0.30 to 0.50 at 0.01 intervals; (c) PAR is simulated with k values that were calculated from the empirical equation[60]; and (d) PAR is simulated with the k values that were calculated from the empirical equation[61].

Figure 7 .
Figure 7.Comparison of the observed and simulated photosynthetically active radiation (PAR) within a crown: (a) PAR is simulated with the k model that was constructed in this study; (b) PAR is simulated with k values ranging from 0.30 to 0.50 at 0.01 intervals; (c) PAR is simulated with k values that were calculated from the empirical equation[60]; and (d) PAR is simulated with the k values that were calculated from the empirical equation[61].

Forests 2018, 9 , 20 Figure 8 .
Figure 8.The daily net photosynthetic production of the whole crown (NPPD) of each pseudowhorl and the whole crown for tree I.The k values for calculating NPPD were respectively determined based on the k model which was constructed in this study (I) and an optimum constant value of 0.33 (II).

Figure 8 .
Figure 8.The daily net photosynthetic production of the whole crown (NPP D ) of each pseudowhorl and the whole crown for tree I.The k values for calculating NPP D were respectively determined based on the k model which was constructed in this study (I) and an optimum constant value of 0.33 (II).

Table 1 .
Summary of parameters of all the photosynthetic light-response (PLR) curves, leaf traits and environmental conditions in 2016-2017.The following leaf variables are given: the maximum net photosynthetic rate (A max ), apparent quantum yield (AQY), dark respiration rate (R d ), photosynthetic light response curve (PLR), relative depth into crown (RDINC), leaf mass per area (LMA) leaf temperature (T leaf ) and vapor pressure deficit (VPD).Mean, Std., Max. and Min.represent the mean value, standard deviation, maximum value and minimum value respectively.The values are based on 752 individual needle-clusters from 36 pseudowhorls belonging to five planted Larix olgensis Henry trees.

Table 3 .
Parameter estimates for the crown photosynthetic light-response (PLR) curves model of planted Larix olgensis trees and the goodness of fit statistics.p-value: the level of significance for the parameter estimations, Std.error: standard error for the parameter estimations, R 2 : coefficient of determination, RMSE: root mean square error, ME: mean error, P: precision estimation.

Table 4 .
Parameter estimates for the extinction coefficient (k) model of planted Larix olgensis trees and the goodness of fit and validation statistics.p-value: the level of significance for the parameter estimations, Std.error: standard error for the parameter estimations, Tree no.: tree number; R 2 : coefficient of determination, RMSE: root mean square error, ME: mean error, P%: precision estimation.