Evaluation of Di ﬀ erent Methods for Estimating the Fraction of Sunlit Leaves and Its Contribution for Photochemical Reﬂectance Index Utilization in a Coniferous Forest

: Proper determinations of light use e ﬃ ciency (LUE) and absorbed photosynthetically active radiation (APAR) are essential for LUE models to simulate gross primary productivity (GPP). This study intended to apply the photochemical reﬂectance index (PRI) to track LUE or APAR variations in a subtropical coniferous forest using tower-based PRI and GPP measurements. To improve the ability of using PRI to track LUE or APAR, a two-leaf approach di ﬀ erentiating sunlit and shaded leaves was used to process the remote sensing and ﬂux data. However, penumbra region, the ‘grey region’ between sunlit and shaded leaves, increases the di ﬃ culty for quantifying the fractions of sunlit and shaded leaves. Firstly, three methods with di ﬀ erent ways on treating the penumbra region were investigated for estimating the fraction of sunlit leaves ( P T ). After evaluating the correlations between observed PRI (PRI obs ) and inversely retrieved PRI (PRI inv ) from estimated P T using the three methods, we found that treating a substantial portion of penumbra region as sunlit leaves was reasonable and using the ratio of canopy reﬂectance to leaf reﬂectance as P T was accurate and e ﬃ cient. Based on this, we used the two-leaf approach to estimate the canopy-level PRI, aiming to evaluate the ability of using PRI as a proxy for LUE or APAR. Results showed that PRI was able to capture half-hourly and daily changes in LUE and APAR, and the two-leaf approach could enhance the correlations between PRI and both LUE and APAR at both half-hourly and daily time steps. Strong diurnal correlations (averaged R = 0.82 from 173 days) between two-leaf PRI and APAR were found on more than 80% days and the relationship between them over the whole study period was also very signiﬁcant ( R 2 > 0.5, p < 0.0001) regardless of di ﬀ erent climate conditions, suggesting that the two-leaf PRI was probably a better proxy for APAR than for LUE at short-term scale as PRI mainly represented the absorbed energy allocated to photoprotection at short time scale and was a direct outcome driven by APAR. However, the scattered relationships of PRI with LUE and APAR indicated there were still many limitations in usage of PRI to accurately estimate physiological parameters a ﬀ ected by changing weather conditions, pigment pool size, etc., which needed further exploration.


Introduction
Remote sensing data have been widely used to calculate gross primary productivity (GPP) in combination with light use efficiency (LUE) models, which commonly express GPP as the product of the amount of absorbed photosynthetically active radiation (APAR) and a LUE term (GPP = LUE × APAR) [1,2]. LUE models are applied to estimation of GPP at the regional or global scale. However, the estimations of the two key input parameters of LUE models, i.e., LUE and APAR, remain challenging [3][4][5]. Proper determination of LUE and APAR are critical for calculating GPP using remote sensing data and LUE models.
In LUE models, stress-based approaches considering a maximum LUE downregulated by environmental factors, such as soil water, vapor pressure deficit (VPD), minimum and maximum air temperature [3,[6][7][8][9][10][11][12]. However, such an assumption probably induces considerable uncertainties in determining GPP [13][14][15][16][17]. As canopy LUE is hard to be directly measured, flux-tower-based canopy LUE is treated as almost the 'real' LUE and is often used for evaluating or validating LUE models at stand and landscape scales [8,[18][19][20][21][22]. It is estimated as GPP derived from net ecosystem CO 2 exchange (NEE) measurements from eddy covariance technique divided by empirically estimated or observed APAR. Previous studies indicated that the fractions of diffuse and direct radiation absorbed by leaves were different [23][24][25], and that LUE was affected by both the quantity and composition of the incoming solar radiation [14,16,26]. With a given value of total incoming radiation, the canopy-level LUE generally increases with increasing fraction of diffuse radiation due to an increase in the canopy fraction that is receiving illumination without photo-saturation [15,17,[26][27][28]. The calculation of APAR is suitable for clear sky conditions due to its theoretic basis on direct radiative transfer regime, but may overestimate APAR of sunlit leaves and underestimate that of shaded leaves when the fraction of diffuse radiation is high.
A promising method used for accurately estimating LUE is the photochemical reflectance index (PRI = (R 531 − R 570 )/(R 531 + R 570 ), where R 531 and R 570 are canopy reflectance at wavelengths 531 and 570 nm, respectively) [29]. The xanthophyll cycle, which is part of the leaf photoprotective mechanisms and functionally related to photosynthetic rates, is closely correlated with the rate of dissipation of excess absorbed energy as heat [30]. During this cycle, the xanthophyll pigment diepoxide violaxanthin is converted rapidly via the intermediate monoepoxide antheraxanthin into the epoxide-free form zeaxanthin by the enzyme violaxanthin de-epoxidase, known as a de-epoxidation sequence [30][31][32][33][34]. As the reflectance change at 531 nm is caused by the appearance of zeaxanthin, PRI is currently considered as the most representative index of the xanthophyll cycle [29,35,36]. As we know, energy absorbed by plants (i.e., APAR) is distributed to three parts (i.e., photochemistry, fluorescence and photoprotection dissipating as heat), and these three portions are competitive. PRI being the proxy for the portion of photoprotection shows the potential to indicate LUE which is an outcome of the competition, as well as APAR which is an input of the competition [22].
Strong correlations were observed between PRI and LUE for individual species, and significant but much weaker correlations were observed for several species as a group [37][38][39][40][41][42][43][44]. However, the relationship between PRI and LUE remains elusive because of the complex influence of canopy architecture [38,42,[45][46][47][48]. Even worse, PRI may be decoupled with LUE when different inherent mechanisms driving PRI and LUE variations under some circumstances [49][50][51]. Under severe stress conditions, some other photoprotection mechanism like photorespiration is also involved in as a sink of excess energy as heat dissipation is not enough for stress protection [52,53]. Decoupling between PRI and photoprotection might occur in response to different sources of stress, mostly reported as severe drought stress and high light intense [49,51], that further make the relationship between PRI and LUE break down [50], limiting the ability of PRI to indicate LUE.
Estimation of APAR of a vegetation canopy is also challenging. Considerable research has been conducted using vegetation indices to calculate the fraction of photosynthetically active radiation absorbed by a canopy (FPAR) and then APAR with photosynthetically active radiation (PAR) mostly from reanalysis remote sensing data [7,54,55]. With products from a moderate resolution imaging spectroradiometer (MODIS), Zhang et al. [56][57][58] presented a method to retrieve the fraction of PAR absorbed by chlorophyll contained throughout the canopy (fPAR chl ) with the advanced radiative transfer model PROSAIL2 and a type of Markov Chain Monte Carlo (MCMC) estimation procedure, which opens a new perspective for accurate estimation of APAR used for photosynthesis. However, estimation of FPAR was sensitive to two structural parameter (i.e., leaf area index (LAI) and the average leaf inclination angle (ALA)) and illumination conditions [59], and estimation of PAR was also of uncertainty [60,61], both of which introduced errors in estimation of APAR. To our limited knowledge, no works have been conducted to estimate APAR directly using remote sensing data. As variation of PRI is directly driven by absorbed energy while LUE changes from more processes [62,63], to some extent PRI is assumed to be a better proxy for APAR than for LUE if the effects of non-physiological factors on PRI are reduced.
To reduce the effects of non-physiological factors on PRI, such as sun-target-view geometry, background reflectance, diffuse sky radiation, and canopy structure, many works have been done in last several years. Hall et al. [64,65] demonstrated that directional changes in observed canopy PRI at a given time interval can be attributed almost entirely to shadow fraction (α s ) variations. The partial derivative of PRI with respect to α s (∂PRI/∂α s ) was used to indicate the relative light use efficiency ∆ε (relative to unstressed, where the unstressed value ∆ε = 1.0, so that ∆ε ranges from 0 to 1) [64,65]. However, this assessment of α s could only be used under clear sky conditions, since the model does not account for diffuse radiation. [66]. The ratio of direct to total irradiance should certainly be considered [26], and simpler and more universal methods by separating the canopy into sunlit and shaded parts under various weather conditions should be developed for LUE models.
Separating sunlit and shaded leaves is vital to accurately predict canopy photosynthesis due to their different light responses [67,68]. Sunlit leaves are illuminated by both the direct and diffuse radiation while shaded leaves are illuminated by only diffuse radiation [69], however, quantifying the fractions of these two parts is hard to achieve. Previous studies could only estimate sunlit and shaded leaf areas of the whole canopy using Chen's method based on solar zenith angle and measurements of leaf area index [23]. Geometric-optical models can be used to simulate the fractions of sunlit and shaded leaves [70][71][72] but need many accurate input parameters and some assumptions. Zhang et al. [73] developed a simple two-leaf approach to separate a canopy-level PRI observation into sunlit and shaded PRI values using a ratio of observed canopy reflectance to leaf reflectance. The two-leaf canopy PRI (PRI t ) was calculated as the sum of sunlit and shaded PRI weighted by their respective sunlit and shaded LAI, based on the principle of linear decomposition of canopy reflectance [70]. Tower-based flux and multi-angle PRI measurements were made over a subtropical forest to develop and validate the algorithm. At both half-hourly and daily time steps, PRI t can effectively improve (>60%) its correlations with LUE derived from the tower flux measurements over the big-leaf PRI taken as the arithmetic average of the multi-angle measurements in a given time interval.
Even though the definitions of sunlit and shaded leaves are clear, except for these two parts, the penumbra region (part sun visible) is also an important element of the canopy since the sun is a solar disc with a finite size instead of a beam with parallel rays [74]. The penumbra region is the 'grey region' between sunlit and shaded leaves, and can be treated as partially sunlit and partially shaded blurring the boundary between sunlit and shaded leaf. Penumbra effect would influence the canopy photosynthesis and stomatal conductance to some extent. The method for estimating the fraction of sunlit leaves presented by Zhang et al. [73] fitted the definition of sunlit leaves; however, this method included a substantial portion of the penumbra region into sunlit leaves, which needed to be verified if the fraction of sunlit leaves is overestimated. This study aimed at testing three methods for estimating the fraction of sunlit leaves that were seen by a remote sensing sensor, which were further used to calculate two-leaf canopy PRI. Comparisons were made among the three methods including two models for estimating the fraction of sunlit leaves based on Zhang's method [73] with different hypotheses on the way to treat penumbra region (hypothesis I: a substantial portion of penumbra region is treated as sunlit leaves; and hypothesis II: a substantial portion of penumbra Remote Sens. 2019, 11, 1643 4 of 22 region suffers a much lower level of light stress than sunlit leaves and similar as shaded leaves), and an geometrical-optical model, which was a novel tool to simulate the fractions of four components of the canopy (i.e., sunlit and shaded leaves and sunlit and shaded backgrounds) [72]. These models were tested against the tower-based flux and PRI datasets in order to find a simple and reliable two-leaf model for improving the ability of PRI measurements for tracking LUE as well as APAR of plant canopies.

Field Measurements
We observed the CO 2 flux and multi-angle canopy spectral reflectance of a sub-tropic mature evergreen coniferous forest stand on a tower at Qianyanzhou Experimental Station (QYZ) in the southeast of China. The CO 2 flux was measured using eddy covariance (EC) method at 39.6 m height flux tower [75]. Half-hourly flux data was provided and processed to derive GPP from measured NEE with daytime ecosystem respiration (R e ) estimated with an empirical equation fitted using nighttime NEE and temperature measurements [76,77], as follows: Multi-angle spectral data were acquired from iAMSPEC system mounted at 31 m on the same tower for EC measurements [73,78], which could concurrently measure incoming solar irradiance and reflected canopy radiance every 2~3 s all day using a dual-channel spectrometer Unispec-DC (PP Systems, Amesbury, MA, USA). A rotator device PTU-D46 (FLIR Systems, Goleta, CA, USA) was used to change view angles to obtain multi-angle canopy reflectance. The view azimuth angle changed at 10 • angular step width with in a range from 45 • to 325 • (defined from geodetic north), and an observation cycle was completed in 15 min with four view zenith angles (firstly fixed at the instantaneous solar zenith angle to finish a horizontal cycle and then at (37 • , 47 • , 57 • ) or (42 • , 52 • , 62 • ) for each view azimuth angle (backward-and forward-looking, changing alternatively) alternatively between two 15-min cycles). After dark current (DC) and calibration of sensors, canopy reflectance (R can ) was calculated as observed canopy reflected radiance divided by solar radiance: where R i and R r are the irradiance and the radiance of the canopy sensor, respectively; the single quote mark means data measured from the reference target. Then PRI was calculated as [73,78]: The canopy-level big-leaf PRI (PRI b ) was simply the arithmetic average of all 2-3 s intervals of observed multi-angle PRI values when the view zenith angle is less than 63 • within half an hour to match the temporal scale of half-hourly LUE.
Available bioclimatic measurements were also made around the same flux tower, mainly including photosynthetically active radiation (PAR, MJ/m 2 /hh), air temperature (Ta, • C), vapor pressure deficit (VPD, kPa), relative humidity (Rh), wind direction, wind speed (m/s) and precipitation (mm). The green leaf area index (LAI) of the canopy was 5.8 for the entire study period, and the clumping index (Ω) was 0.57 estimated from hemispherical photos for this forest. Data from day of year (DOY) 101 to 275 in 2013 were used in this study.

Calculations of APAR and LUE
APAR was calculated using the algorithm in a two-leaf LUE model (TL-LUE) developed by He et al. [15]. This algorithm differentiates transfers of diffuse and direct radiation within a canopy. Sunlit leaves absorb both direct and diffuse radiation, while shaded leaves are only exposed to diffuse radiation. The total LAI of the canopy can be separated into sunlit LAI (L sun ) and shaded LAI (L sh ) according to half-hourly solar zenith angle (θ) based on the methods presented by Chen et al. [23,24]. APAR on per unit leaf area was calculated separately for sunlit (APAR sun ) and shaded leaves (APAR sh ) (detailed derivation can be found in Appendix A), hence the two-leaf APAR and LUE were calculated as [15,73,78]:

Two-Leaf PRI Model
Treating the canopy as two big leaves, i.e., sunlit and shaded leaves, and considering their distinct difference of absorbed radiation and PRI, a two-leaf approach was developed to derive PRI t using multi-angle observations.
To derive fractions of observed sunlit and shaded leaves, observed canopy reflectance (R can ) for each view angle was separated to four components on the basis of a four-scale bidirectional reflectance model [70], as following: where P T and P S are the fractions of sunlit and shaded leaves, respectively; P G and P Z are the fractions of sunlit and shaded background, respectively; R T and R S are the reflectance of sunlit and shaded leaves, respectively; R G and R Z are the reflectance of sunlit and shaded background, respectively. The fractions of four components are irrelevant with respect to wavelength, whereas reflectance of these four elements is wavelength dependent.

Estimations of Observed Fraction of Background
The observed fraction of background (P VG ) at a given view angle (θ v ) is composed of P G and P Z , which are hard to be estimated separately from observations. Assuming the background had no vegetation and its reflectance was constant, P VG , which was also the gap fraction from the viewer direction, could be totally calculated based on the methods presented by Chen et al. [70] as: 2.3.2. Estimations of P T and P S According to Equation (8), the sum of the fractions of sunlit and shaded leaves and sunlit and shaded backgrounds is unit, thus the fraction of shaded leaves (P S ) was: Estimation of canopy P T was based on both canopy and leaf reflectances presented by Zhang et al. [73]. The ratio of the observed canopy reflectance (R can ) to leaf reflectance (R leaf ) at 670 nm (the multi-scattering and the contribution from background in this band is the lowest), defined as R SL , was taken to be approximately equal to the degree to which the canopy was exposed to radiation. Determination of the observed P T at each view angle was tested on two different strategies of whether or not considering the contribution to canopy reflectance from shaded leaves, which were only exposed to diffuse radiation. A hybrid geometric-optical model (GOST) with ray tracing method was also tested to estimate P T , which was also treated as a reference to exam P T calculated from the other two methods, as GOST was a novel model to simulate the fractions of the four components of a canopy presented in equation (8).
Model I: Treating P T as the fraction of the canopy exposed to total radiation.
According to the definition of sunlit leaf, a substantial portion of penumbra region was treated as sunlit leaves. Hypothesizing this part of penumbra region suffered the same level of light stress on photosynthesis as sunlit leaves, without separating direct and diffuse radiation, P T in this model (P T,I ) was equal to R SL : This model could represent the total degree of leaves exposed to radiation. It satisfied the assumption that photosynthesis often saturates under high radiation while it was supposed to be non-saturate under diffuse radiation. However, as the canopy reflectance was composed with reflectances of both sunlit and shaded leaves, R SL might include the contribution from shaded leaves of a part of penumbra region. Therefore, this model was hypothesized overestimating the real P T .
Model II: Treating P T as the fraction of the canopy exposed to total radiation subtracted by the fraction of shaded leaves exposed to diffuse radiation.
Hypothesizing a substantial portion of penumbra region suffered a much lower level of light stress than sunlit leaves and similar to shaded leaves, we further hypothesized that the direct radiation contributes to the light stress while the diffuse radiation cause none light stress considering the penumbra effect. Therefore, P T in this model (P T,II ) was the degree of the canopy exposed to direct radiation, and could be treated as the ratio of canopy reflectance from direct incident radiation (instead of total incident radiation) to leaf reflectance. With the ratio of diffuse radiation to total radiation (R dif ), the fraction of shaded leaves was subtracted from R SL to calculate P T,II . The canopy reflectance come from direct and diffuse radiation, can be expressed as: where PAR dif is diffuse PAR, of which calculation can be found in Appendix A; R can,dir and R can,dif are canopy reflectances from direct and diffuse incident radiation, respectively. Therefore, As R can can be divided into R can,dir and R can,dif by the direct and diffuse ratios to total radiation, R can,dif was: thus, with equation (10), Finally, P T,II can be finally expressed as This model also satisfied the definition that sunlit leaves absorbed both direct and diffuse radiation while shaded leaves only absorbed diffuse radiation. However, as the real sunlit and shaded leaves were not always easy to distinguish, this model treated a substantial portion of penumbra region as shaded leaves, hence we hypothesized that it may over-calculate the fraction of shaded leaves and hence underestimate the real P T .
Model III: Modeling P T using a hybrid geometric-optical model with ray tracing method for sloping terrains (GOST) [71,72,79].
In the geometric-optical model for sloping terrains (GOST), the shape of trees was simplified as adiactinic cylinder and cone, and there was no penumbra region within the canopy. It assumed a Neyman distribution, also called double-Poisson distribution, of the trees. A simplified ray tracing method was designed for simulating the ratio of the viewed sunlit leaves to all viewed leaves (S PT ). The basic idea of this method was first to describe the foliage spatial and angular distributions, and then to penetrate a view line into the canopy. If the view line does not touch any leaves in the forest scene, it was not considered. Otherwise, the first intersection point (FIP) of the view line and leaves could be found. If there was no other leaves shaded the FIP in the sunlight direction, it was a sunlit point in the view direction. If there were other leaves between the sun's position and the FIP, the FIP was the shaded point. Through repeating these steps many times, the probability of viewing sunlit leaves could be separated from the viewed leaves.
After all the ray tracing procedures, the percentage of sunlit points S PT reached by the view lines in the view direction could be simulated. The fraction of sunlit leaves in the canopy was To drive GOST, several canopy structure parameters were measured and presented in Table 1. Table 1. The canopy structure parameters for the geometric-optical model for sloping terrains (GOST).
Ha is the height of the lower part of the tree (trunk space), Hb is the height of crown, r is the diameter of the crown, α is the half apex angle, Ws is the mean width of elements shadows cast inside tree crowns, γ E is the needle-to-shoot area ratio and G(θ) is the projection of unit leaf areas.

Estimation of Two-Leaf PRI
As PRI of the background with no vegetation was presumed to be constant [80] and the light condition was supposed to be steady, radiation and PRI of sunlit (PRI sun ) and shaded (PRI sh ) leaves were assumed to be constant within 15 min. Thus, PRI obs at different view angles varies only because of the variable sunlit and shaded portions of the canopy, and then could be separated into two components: PRI sun and PRI sh were estimated with multi-angle observations acquired within 15 min using least squares regression.
Finally, PRI t was calculated with the weights of L sun and L sh :

Statistical Data Analysis
All data measured with the solar zenith angle less than 75 • and the view zenith angle less than 63 • were selected to reduce the influence of low quality of measurements when solar and view zenith angle were too large. In total, there were nearly one hundred thousand of 2-3 s intervals observed samples of multi-angle reflectance and PRI obtained under different weather conditions. To test the three methods for estimating the fraction of sunlit leaves, samples within each 15-min observation were separately analyzed. Once we derived PRI sun , PRI sh , P T and P S from the three models, PRI values at each view angle were inversely calculated (PRI inv ) using reversed Equation (19), and then compared with each observed PRI (PRI obs ) within 15 min. The linear correlations coefficients between PRI inv and PRI obs within 15-min intervals were calculated. The frequencies of correlations coefficients at different significant levels were used to indicate the accuracy of the three models.
Correlations between PRI and LUE or APAR at both half-hourly and daily temporal scales were presented to evaluate the ability of PRI t to indicate LUE and APAR against that of PRI b . For diurnal relationship analysis, half-hourly PRI, LUE and APAR obtained from 7:00 to 17:30 local time were selected to process linear regression. Half-hourly PRI was averaged from two 15-min PRI. For the whole study period, the linear and non-linear (exponential) relationships of PRI with LUE and APAR were assessed for discussing the utilization of PRI to indicate LUE and APAR at both half-hourly and daily scales. Statistical analyses were performed using MATLAB and Microsoft Excel software.

Test of Three Methods for Estimating the Fraction of Sunlit Leaves
With the three models, three sets of P T were derived from multi-angle datasets. Figure 1 showed the variations of three sets of P T with different view angles within an observation cycle from 10:45 to 11:00 on April 11. P T from model I and model II showed quite similar trends but totally different values. Considering the penumbra effect, a substantial portion of penumbra was considered as shaded leaves in model II, leading to the much lower values of P T,II than that of P T,I . P T from model III exhibited a larger variation range than the other two models, with a slightly different variation trend. High values of P T,III matched those of P T,I well, while low values were close to those of P T,II . High values of P T were observed from backward looking (sensor pointing north) and low values were observed from forward looking (sensor pointing south). At the same view azimuth angle, larger view zenith angle had bigger values of P T . P T,III from forward looking had low values with changing view zenith angle. Model III assumed a uniform canopy and used statistical methods to simulate P T [70,79], which meant that it needed a large sample for satisfying this assumption. This indicated that model III cannot capture small structural changes from the real canopy, which was heterogeneous, so that the variation of P T,III was different from the other two. accuracy of the three models.
Correlations between PRI and LUE or APAR at both half-hourly and daily temporal scales were presented to evaluate the ability of PRIt to indicate LUE and APAR against that of PRIb. For diurnal relationship analysis, half-hourly PRI, LUE and APAR obtained from 7:00 to 17:30 local time were selected to process linear regression. Half-hourly PRI was averaged from two 15-minute PRI. For the whole study period, the linear and non-linear (exponential) relationships of PRI with LUE and APAR were assessed for discussing the utilization of PRI to indicate LUE and APAR at both half-hourly and daily scales. Statistical analyses were performed using MATLAB and Microsoft Excel software.

Test of Three Methods for Estimating the Fraction of Sunlit Leaves
With the three models, three sets of PT were derived from multi-angle datasets. Figure. 1 showed the variations of three sets of PT with different view angles within an observation cycle from 10:45 to 11:00 on April 11. PT from model I and model II showed quite similar trends but totally different values. Considering the penumbra effect, a substantial portion of penumbra was considered as shaded leaves in model II, leading to the much lower values of PT,II than that of PT,I. PT from model III exhibited a larger variation range than the other two models, with a slightly different variation trend. High values of PT,III matched those of PT,I well, while low values were close to those of PT,II. High values of PT were observed from backward looking (sensor pointing north) and low values were observed from forward looking (sensor pointing south). At the same view azimuth angle, larger view zenith angle had bigger values of PT. PT,III from forward looking had low values with changing view zenith angle. Model III assumed a uniform canopy and used statistical methods to simulate PT [70,79], which meant that it needed a large sample for satisfying this assumption. This indicated that model III cannot capture small structural changes from the real canopy, which was heterogeneous, so that the variation of PT,III was different from the other two. The correlations between observed PRI (PRIobs) and inversely retrieved PRI (PRIinv) from estimated PT were used to evaluate the accuracy of three models. The frequency of the correlations at different significance levels was a proxy for accuracy. Figure. 2 showed variations of the The correlations between observed PRI (PRI obs ) and inversely retrieved PRI (PRI inv ) from estimated P T were used to evaluate the accuracy of three models. The frequency of the correlations at different significance levels was a proxy for accuracy. Figure 2 showed variations of the frequencies of the models I and II over the study period (number of measurements = 5700). PRI inv from both models matched observed PRI well (>70% samples) with a significance level of 95% over the entire period. Model I performed marginally better than model II either in each month or over the entire period. The solid lines in Figure 2 showed variations of the mean residuals (multiplied by 100) throughout the six months, which include estimation errors and PRI of sunlit background (shaded background was approximately 0), as PRI of background was presumed to be constant and not taken into account. Residuals in model I were considerably smaller than those in model II, indicating model I was better than model II. Moreover, the difference of PRI between sunlit and shaded background in the dry season (from July to September) [78] was theoretically higher than that of wet background in the rainy season (from April to June), indicating that residuals should be large in the dry season [73]. In model I, mean residuals in July and August were larger than other months. The mean residual in September was small probably due to large solar zenith angle and more litter cover on the ground. However, in model II, a big part of the penumbra region was treated as shaded leaves, so that the fraction of leaves exposed to diffuse radiation contributed to PS. Residuals in model II varied in a large range and were much larger in the rainy season than in the dry season, proving that model II was of lower accuracy due to underestimated PT especially in the rainy season when the fraction of diffuse radiation was high. As in model II, fThis further indicated that the penumbra region was better represented by sunlit leaf conditions under light stress, and model I had a higher accuracy.
As the ray tracing treatment of GOST was computationally demanding, model III was more computationally intensive than the other two models by several orders of magnitude to simulate PT and PS. Therefore, only data acquired in April were processed using model III to compare its accuracy with those of the other two models. Figure 3 showed differences between the frequencies of the three models in April. The accuracy of model III was higher than the other two models under most significant levels. However, its leading position faded away with increases of significance level from 0.05 to 0.001. Its accuracy was lower than that of model I at significance level of 0.0001. The mean residuals (multiplied by 1000) of PRIinv of model III were larger than residuals of model I but smaller than those of model II. The solid lines in Figure 2 showed variations of the mean residuals (multiplied by 100) throughout the six months, which include estimation errors and PRI of sunlit background (shaded background was approximately 0), as PRI of background was presumed to be constant and not taken into account. Residuals in model I were considerably smaller than those in model II, indicating model I was better than model II. Moreover, the difference of PRI between sunlit and shaded background in the dry season (from July to September) [78] was theoretically higher than that of wet background in the rainy season (from April to June), indicating that residuals should be large in the dry season [73]. In model I, mean residuals in July and August were larger than other months. The mean residual in September was small probably due to large solar zenith angle and more litter cover on the ground. However, in model II, a big part of the penumbra region was treated as shaded leaves, so that the fraction of leaves exposed to diffuse radiation contributed to P S . Residuals in model II varied in a large range and were much larger in the rainy season than in the dry season, proving that model II was of lower accuracy due to underestimated P T especially in the rainy season when the fraction of diffuse radiation was high. As in model II, fThis further indicated that the penumbra region was better represented by sunlit leaf conditions under light stress, and model I had a higher accuracy.
As the ray tracing treatment of GOST was computationally demanding, model III was more computationally intensive than the other two models by several orders of magnitude to simulate P T and P S . Therefore, only data acquired in April were processed using model III to compare its accuracy with those of the other two models. Figure 3 showed differences between the frequencies of the three models in April. The accuracy of model III was higher than the other two models under most significant levels. However, its leading position faded away with increases of significance level from 0.05 to 0.001. Its accuracy was lower than that of model I at significance level of 0.0001. The mean residuals (multiplied by 1000) of PRI inv of model III were larger than residuals of model I but smaller than those of model II. and PS. Therefore, only data acquired in April were processed using model III to compare its accuracy with those of the other two models. Figure 3 showed differences between the frequencies of the three models in April. The accuracy of model III was higher than the other two models under most significant levels. However, its leading position faded away with increases of significance level from 0.05 to 0.001. Its accuracy was lower than that of model I at significance level of 0.0001. The mean residuals (multiplied by 1000) of PRIinv of model III were larger than residuals of model I but smaller than those of model II. The enumbra region occupied a substantial portion of the conifer forest canopy. To accurately estimate the fraction of sunlit leaves in the canopy, the way to treat penumbra region was vital for the estimation accuracy [74]. Models I and II treated penumbra region differently, resulting in quite different P T values ( Figure 1). Model III made some assumptions so that penumbra region does not exist in the model. Since there was no useful method to quantify the penumbra region for now, we made two different hypotheses for two models to justify how we should make the boundary between sunlit and shaded leaves clear. By testing the accuracy of inversely retrieved PRI from estimated using three models, we found that model II did not perform as well as the other two models. Model III was generally the best one of higher accuracy compared to model I, which indicated that GOST was a novel tool to simulate the fractions of the four components of the canopy. However, model III treated the canopy as a uniform surface which caused some differences from the real trees distribution in simulating P T and P S at different view angles. Furthermore, model III was very computationally demanding and required six months to process six months of data using a high grade workstation while model I was computationally efficient only needing tens of minutes with comparable high accuracy. On the other hand, even though model I may introduce some errors from hypotheses [73], the similar accuracy of model I with model III proved the effectiveness of the method to retrieve two-leaf canopy PRI [73]. Results presented by Middleton et al. [81] indicated the PRI values of penumbra region (sunlit-shaded mixed group) were more closer to sunlit leaves than shaded leaves. Likewise, our results also demonstrated that the penumbra region mostly should be treated as sunlit leaves. Therefore, we found model I is the optimal two-leaf approach to calculate canopy-level PRI.

Evaluations of the Two-Leaf PRI (PRIt) to Indicate LUE and APAR
Firstly, PRI t retrieved using model I was analyzed with LUE and APAR, and compared with PRI b . Figure 4 displayed the diurnal variations of half-houly LUE, APAR, PRI b and PRI t averaged from the whole study period and the coefficients of determination (R 2 ) of LUE and APAR with PRI b and PRI t . The grey areas denote the standard deviation of each of the four terms. Generally, LUE and PRI t showed similar variation trends, while APAR had an opposite changing trend. PRI b showed a quite different diurnal variation that increased in early morning and changed marginally during the rest of day. PRI t could well track changes of both LUE and APAR, especially APAR. PRI b could also track change of LUE, but hardly be able to track change of APAR. As observed PRI b at early morning and late afternoon mostly represented leaves of upper canopy (mostly sunlit leaves), of which values were relatively lower than those of lower canopy [73], PRI b was even lower than that at midday owing to sun-viewer geometry. The standard deviation of PRI b was larger than that of PRI t , indicating that PRI b was affected by non-physiological factors, especially the sun-viewer geometry, and PRI t could reduce these effects [73]. In early morning and late afternoon, a large part of sunlit leaves on top of the canopy, which was of low PRI values, was observed due to the sun-viewer geometry, so that the value of PRI b was as low as that around noon and had marginal diurnal variation. PRI t reduced effect of sun-viewer geometry and had similar diurnal variation with LUE. Seasonal variations of daily PRIb, PRIt, LUE and APAR averaged from half-hourly data from 10:00 to 15:30 per day were shown in Figure 5. LUE and APAR exhibited quite inversely similar seasonal trends. During the period from day 218 to 225, when a long-term heat wave with high values of APAR occurred [78], both PRIb, PRIt and LUE decreased, whereas PRIb seriously declined but the relationship between PRIt and LUE also broken down under severe heat wave [73]. As the study period was the growing season of this conifer forest while the chlorophyll content and LAI changed marginally [73], there was no obvious seasonal variations of these four terms were observed.   Seasonal variations of daily PRI b , PRI t , LUE and APAR averaged from half-hourly data from 10:00 to 15:30 per day were shown in Figure 5. LUE and APAR exhibited quite inversely similar seasonal trends. During the period from day 218 to 225, when a long-term heat wave with high values of APAR occurred [78], both PRI b , PRI t and LUE decreased, whereas PRI b seriously declined but the relationship between PRI t and LUE also broken down under severe heat wave [73]. As the study period was the growing season of this conifer forest while the chlorophyll content and LAI changed marginally [73], there was no obvious seasonal variations of these four terms were observed. Seasonal variations of daily PRIb, PRIt, LUE and APAR averaged from half-hourly data from 10:00 to 15:30 per day were shown in Figure 5. LUE and APAR exhibited quite inversely similar seasonal trends. During the period from day 218 to 225, when a long-term heat wave with high values of APAR occurred [78], both PRIb, PRIt and LUE decreased, whereas PRIb seriously declined but the relationship between PRIt and LUE also broken down under severe heat wave [73]. As the study period was the growing season of this conifer forest while the chlorophyll content and LAI changed marginally [73], there was no obvious seasonal variations of these four terms were observed.     Figure 6 showed the number of days on which the correlation coefficients (R) were significant at different levels (p). The R between half-hourly PRI b or PRI t and LUE (Figure 4a) or APAR (Figure 6b) on individual days were calculated with more than five good-quality observations from 7:00 to 17:30, when the solar zenith angle was less than 75 • . As PRI was known to be positively correlated with LUE during the growing season, the numbers of days that PRI t was positively correlated to LUE at the significance levels of 0.01 and 0.001 are almost double those of PRI b (Figure 6a). The two-leaf treatment obviously improved the ability of PRI to track diurnal variations of LUE in comparison to the big-leaf PRI. The correlation between PRI t and APAR (Figure 6b) was highly enhanced as among 173 days, with more than 80% days that PRI t was highly significantly correlated to APAR (p < 0.001) while more than 60% days that PRI b was weakly related to APAR (p > 0.05). These results suggested that two-leaf PRI was a better proxy for APAR at a short temporal scale than LUE.
Remote Sens. 2018, 10, x FOR PEER REVIEW  12 of 22 17:30, when the solar zenith angle was less than 75°. As PRI was known to be positively correlated with LUE during the growing season, the numbers of days that PRIt was positively correlated to LUE at the significance levels of 0.01 and 0.001 are almost double those of PRIb (Figure. 6a). The two-leaf treatment obviously improved the ability of PRI to track diurnal variations of LUE in comparison to the big-leaf PRI. The correlation between PRIt and APAR ( Figure. 6b) was highly enhanced as among 173 days, with more than 80% days that PRIt was highly significantly correlated to APAR (p<0.001) while more than 60% days that PRIb was weakly related to APAR (p>0.05). These results suggested that two-leaf PRI was a better proxy for APAR at a short temporal scale than LUE.  Figure 7 illustrated the relationships of PRIb and PRIt to APAR and LUE, respectively, at half-hourly temporal scales with observations from 7:00 to 17:30 throughout the study period in all weather conditions. PRIb was less correlated with either APAR or LUE than PRIt. The relationship between PRIt and APAR was much stronger than that between PRIb and APAR with significantly increased R 2 value (increased ~71%; Figure 7a, 7c). Generally, these results again proved the feasibility of the two-leaf approach tested in this study and presented by Zhang et al. [73]. In addition, both PRIb and PRIt showed higher correlation with APAR than LUE, which indicated that PRI was probably a better proxy for APAR than for LUE. Even though previous studies were focused on using PRI to indicate LUE, stronger correlation between PRI and APAR than that between PRI and LUE at short-term across a growing season have been reported for different forests by Soudani et al. [44]. Decrease of PRI with increased APAR is linked with an increase in light absorption associated with the conversion of violaxanthin into antheraxanthin and zeaxanthin pigments, which can be detected by the decrease in reflectance at 531 nm [35]. PRI is more sensitive to absorbed energy, exhibiting closely increasing with APAR even when net photosynthesis rate and LUE are still low (often occur with low stomatal conductance in the afternoon) [49]. APAR are significant at different levels (p). The value of R between half-hourly PRI b /PRI t and LUE and APAR was calculated per day with more than five good-quality observations from 7:00 to 17:30. Figure 7 illustrated the relationships of PRI b and PRI t to APAR and LUE, respectively, at half-hourly temporal scales with observations from 7:00 to 17:30 throughout the study period in all weather conditions. PRI b was less correlated with either APAR or LUE than PRI t . The relationship between PRI t and APAR was much stronger than that between PRI b and APAR with significantly increased R 2 value (increased~71%; Figure 7a,c). Generally, these results again proved the feasibility of the two-leaf approach tested in this study and presented by Zhang et al. [73]. In addition, both PRI b and PRI t showed higher correlation with APAR than LUE, which indicated that PRI was probably a better proxy for APAR than for LUE. Even though previous studies were focused on using PRI to indicate LUE, stronger correlation between PRI and APAR than that between PRI and LUE at short-term across a growing season have been reported for different forests by Soudani et al. [44]. Decrease of PRI with increased APAR is linked with an increase in light absorption associated with the conversion of violaxanthin into antheraxanthin and zeaxanthin pigments, which can be detected by the decrease in reflectance at 531 nm [35]. PRI is more sensitive to absorbed energy, exhibiting closely increasing with APAR even when net photosynthesis rate and LUE are still low (often occur with low stomatal conductance in the afternoon) [49].
The relationships of daily mean PRI b and PRI t to APAR and LUE averaged from observations between 10:00 and 15:30 throughout the study period were shown in Figure 8. Significant linear correlations were also found for all relationships, while the relationships of PRI b and PRI t with LUE were more non-linear (exponential). The correlation between PRI t and LUE was about 35% stronger than that between PRI b and LUE, although some low PRI t values corresponding to high LUE (Figure 8b,d). The relationship between daily PRI t and APAR was also stronger than that between daily PRI b and APAR (Figure 8a,c), but was weaker than the relationship at half-hourly scale (Figure 7c). The relationships of daily mean PRIb and PRIt to APAR and LUE averaged from observations between 10:00 and 15:30 throughout the study period were shown in Figure 8. Significant linear correlations were also found for all relationships, while the relationships of PRIb and PRIt with LUE were more non-linear (exponential). The correlation between PRIt and LUE was about 35% stronger than that between PRIb and LUE, although some low PRIt values corresponding to high LUE (Figure 8b, 8d). The relationship between daily PRIt and APAR was also stronger than that between daily PRIb and APAR (Figure 8a, 8c), but was weaker than the relationship at half-hourly scale (Figure 7c). Even though the ability of PRI to indicate LUE had been discussed in previous works [4,22,42,44,47,64,78,82], PRI was still difficult to be use to estimate LUE over different tempo-spatial scales as their relationship was influenced by many factors, such as pigment concentration, sun-target-view geometry, background reflectance, diffuse sky radiation, canopy structure. This limited the effectiveness of PRI to be involved in LUE models for GPP estimations [2,41,43,46,65,78,83,84]. Generally, either PRI b or PRI t was more linked to APAR than LUE at different temporal scales (Figures 7 and 8). Porcar-Castell et al. [62] demonstrated that APAR was directly linked with dissipation of a part of energy as heat, which was related to PRI [29,33,36,85], while LUE was dominated by further processes linked with energy distribution. This indicated the stronger mechanistic linkage between PRI and APAR than LUE. PRI t could better capture variations of half-hourly APAR (Figure 7) than those of daily APAR, indicating that the two-leaf PRI was efficient to track short-term (sub-day or even instantaneous) variations of APAR. The penumbra region probably contributed a large amount of absorbed radiation almost like or even more than that absorbed by sunlit leaves in this conifer forest, demonstrating that treating penumbra region as sunlit leaves was suitable for understanding the link between fast-change APAR and canopy PRI (PRI t ). Mechanistically, PRI mainly represented the absorbed energy allocated to photoprotection at short time scale [62,73,81,84,[86][87][88], and was consequently effective to indicate APAR under generally steady conditions (when the distribution of absorbed energy does not change dramatically) as PRI was a direct outcome driven by APAR. But more importantly, this result demonstrated that, at short term, PRI was only able to track changes in LUE driven by APAR, which also had been mentioned by Zhang et al. [73]. This finding further revealed the limitation of using PRI to track diurnal variation of LUE at some circumstances that LUE was not mainly driven by radiation. Even though the ability of PRI to indicate LUE had been discussed in previous works [4,22,42,44,47,64,78,82], PRI was still difficult to be use to estimate LUE over different tempo-spatial scales as their relationship was influenced by many factors, such as pigment concentration, sun-target-view geometry, background reflectance, diffuse sky radiation, canopy structure. This limited the effectiveness of PRI to be involved in LUE models for GPP estimations [2,41,43,46,65,78,83,84]. Generally, either PRIb or PRIt was more linked to APAR than LUE at different temporal scales (Figures 7, 8). Porcar-Castell et al. [62] demonstrated that APAR was directly linked with dissipation of a part of energy as heat, which was related to PRI [29,33,36,85], while LUE was dominated by further processes linked with energy distribution. This indicated the stronger mechanistic linkage between PRI and APAR than LUE. PRIt could better capture variations of half-hourly APAR (Figure 7) than those of daily APAR, indicating that the two-leaf PRI was efficient to track short-term (sub-day or even instantaneous) variations of APAR. The penumbra region probably contributed a large amount of absorbed radiation almost like or even more than that absorbed by sunlit leaves in this conifer forest, demonstrating that treating penumbra region as sunlit leaves was suitable for understanding the link between fast-change APAR and canopy PRI (PRIt). Mechanistically, PRI mainly represented the absorbed energy allocated to photoprotection at short time scale [62,73,81,84,[86][87][88], and was consequently effective to indicate APAR under generally steady conditions (when the distribution of absorbed energy does not change dramatically) as PRI was a direct outcome driven by APAR. But more importantly, this result demonstrated that, at short term, PRI was only able to track changes in LUE driven by APAR, which also had been mentioned by Zhang et al. [73]. This finding further revealed the limitation of using PRI to track diurnal variation of LUE at some circumstances that LUE was not mainly driven by radiation. Figure 8. Relationships of PRI b and PRI t with APAR and LUE at daily temporal scales averaged from half-hourly data from 10:00 to 15:30 throughout the study period. The two R 2 represent the coefficients of determination of linear (red) and exponential (green) regressions, respectively. Above all, the two-leaf approach, which used the canopy reflectance ratio to leaf reflectance to represent the fraction of sunlit leaves, was simple and very effective, verifying the work presented by Zhang et al. [73]. During the process of two-leaf PRI calculation, it was proved that the penumbra region should be mostly treated as sunlit leaves only on concept. If intending to quantify the penumbra effect, light detection and ranging (LiDAR) could be a useful technology to simulate the radiation distribution within the canopy, and then identify the sunlit, penumbra and shaded elements, even could map the 3D distribution of LUE within the canopy [89]. PRI could be used to indicate both of the two key parameters (i.e., LUE and APAR) of LUE model, while two-leaf PRI was a better proxy for APAR at a shorter time scale. It was promising to estimate APAR from hyperspectral data through PRI t with high spatial and temporal resolution, while to estimate LUE was limited under the circumstance that LUE was mainly driven by radiation. However, to generalize the utilization of PRI, further works are needed for clarifying the benefits and limitations of PRI in indicating photosynthetic parameters. To promote using PRI as a proxy for LUE, their relationships at different temporal scales and under different weather conditions should be investigated in different ecosystems.

Conclusions
Shaded leaves, as well as sunlit leaves, contribute quite differently to reflectance and GPP at canopy level. Therefore, separating the canopy into sunlit and shaded parts for estimation of GPP using remote sensing technology is very necessary. Penumbra effect increases the difficulty for quantifying the fractions of sunlit and shaded leaves. In this paper, we investigated three methods with different hypotheses on treatments of the penumbra region for estimating the fraction of sunlit leaves. Inversely retrieved PRI (PRI inv ) from estimated P T were compared to multi-angle observed PRI (PRI obs ) to evaluate the accuracy of the three methods. The high accuracy of model I demonstrated that the penumbra region should be mostly treated as sunlit leaves, of which the fraction could be simply estimated using the ratio of canopy reflectance to leaf reflectance.
The simple two-leaf approach presented by Zhang et al. [73] was used to process the remote sensing and flux data, aiming to evaluate the ability of using PRI as a proxy for LUE or APAR. Generally, PRI was able to capture half-hourly and daily changes in LUE as well as APAR. At either half-hourly or daily time steps, the two-leaf approach significantly enhanced the correlations between PRI and LUE and APAR. The correlation between two-leaf PRI and APAR was very significant under all weather conditions, especially at half-hourly temporal scale. These results suggested that two-leaf PRI was probably a better proxy for APAR than LUE. However, in order to make PRI effective for understanding plant photosynthetic status under different weather conditions, a lot of works need to be done before a general method for reducing the effects of non-physiological factors on PRI can be wildly used.

Patents
Chinese patent: An approach for retrieving canopy photochemical reflectance index of separated sunlit and shaded leaves by using multi-angle observation data. Authorized: ZL201310283751.4

Acknowledgments:
We gratefully acknowledge the support from the Qianyanzhou station of the Chinese Academy of Sciences for providing flux and meteorology datasets. We further acknowledge Jing M. Chen (University of Toronto, Canada) for his valuable discussion on different methods for estimation of the fraction of sunlit leaves. Many thanks to three anonymous reviewers who provided helpful comments to the improvement of the manuscript.

Abbreviations
The following abbreviations are used in this manuscript: The fraction of PAR absorbed by chlorophyll G(θ) The projection of unit leaf areas GPP Gross primary production GOST Geometric-optical model Ha The height of the lower part of the tree (trunk space) Hb The The observed fraction of sunlit leaves P S The observed fraction of shaded leaves P VG The observed fraction of background P G The observed fraction of sunlit background P Z The observed fraction of shade background QYZ Qianyanzhou Experimental Station r The diameter of the crown R 0 The extraterrestrial radiation at the top of the atmosphere R can Canopy reflectance R T Reflectance of sunlit leaves R G Reflectance of sunlit background R S Reflectance of shaded leaves R Z Reflectance of shaded background R dif The ratio of diffuse radiation to total radiation R e Daytime ecosystem respiration R g The global solar radiation on the surface of the earth R leaf Leaf reflectance R SL The ratio of canopy reflectance to leaf reflectance Rh

Relative humidity S PT
The ratio of the sunlit point in a tree Ta Air temperature VPD Vapor pressure deficit W S The mean width of elements shadows cast inside tree crowns

Appendix A
Estimation of PAR dif , PAR dir , APAR sun and APAR sh :APAR on per unit leaf area is calculated separately for sunlit (APAR sun ) and shaded leaves (APAR sh ): where α is the albedo related to vegetation types setting as 0.15 for coniferous forest; PAR dif and PAR dir are the diffuse and direct components of incoming PAR, respectively, and they are calculated using equation (A6); PAR dif,u is the diffuse PAR under the canopy; (PAR dif − PAR dif,u )/LAI represents the diffuse PAR on per unit leaf area within the canopy; C quantifies the contribution of multiple scattering of the total PAR to the diffuse irradiance per unit leaf area within the canopy; β is mean leaf-sun angle and set as 60º for a canopy with a spherical leaf angle distribution; δ is a correction for nonlinear response of leaf photosynthesis to the vertical variation of diffuse radiation within the canopy that makes APAR sh represents absorbed radiation times light response of photosynthesis; and θ is the solar zenith angle. Diffuse and direct PAR are partitioned [15,23,24] with parameters calibrated using daily clearness index (CI) and total incoming radiation as: PAR di f = PAR × 0.7702 + 3.6895CI − 15.4527CI 2 + 16.9828CI 3 − 5.7773CI 4 (A6) CI is a ratio of the global solar radiation on the surface of the earth (R g ) and the extraterrestrial radiation at the top of the atmosphere (R 0 ). R 0 for any day and any moment can be calculated using solar zenith (θ) at a given place as: R 0 = S 0 × (1 + 0.033 cos(360DOY/365)) × cos θ (A7) where S 0 is the solar constant (1367 W/m 2 ); DOY denotes day of the year. With observed half-hourly total incident radiation (R g ), CI is calculated as R g /R 0 .