Impacts of Leaf Age on Canopy Spectral Signature Variation in Evergreen Chinese Fir Forests

Significant gaps exist in our knowledge of the impact of leaf aging on canopy signal variability, which limits our understanding of vegetation status based on remotely sensed data. To understand the effects of leaf aging at the leaf and canopy scales, a combination of field, remote-sensing and physical modeling techniques was adopted to assess the canopy spectral signals of evergreen Cunninghamia forests. We observed an approximately 10% increase in Near-Infrared (NIR) reflectance for new leaves and a 35% increase in NIR transmittance for mature leaves from May to October. When variations in leaf optical properties (LOPs) of only mature leaves, or both new and mature leaves were considered, the Geometric Optical and Radiative Transfer (GORT) model-simulated canopy reflectance trajectory was more consistent with Landsat observations (R2 increased from 0.37 to 0.82~0.89 for NIR reflectance, and from 0.35 to 0.67~0.88 for EVI2, with a small RMSE (0.01 to 0.02)). This study highlights the importance of leaf age on leaf spectral signatures, and provides evidence of age-dependent LOPs that have important impacts on canopy reflectance in the NIR band and EVI2, which are used to monitor canopy dynamics and productivity, with important implications for RS and forest ecosystem ecology.


Introduction
Forests cover approximately 30% of the Earth's land area (4.2 × 10 9 hectares).Globally, forests play critical roles in providing goods and services for terrestrial ecosystems, including filtering water, controlling water runoff, protecting soil, regulating the climate, and cycling and storing nutrients [1,2].Many important biophysical processes in forests are conducted through leaves, including photosynthesis, transpiration, respiration, and light interception.Forests and other land vegetation currently remove approximately 30% of anthropogenic CO 2 emissions from the atmosphere through photosynthesis [3][4][5].While high value has been placed on remote sensing (RS) for ecological research, management and modeling of forest canopy status at an ecosystem scale, a concomitant increase in understanding the factors that affect canopy reflectance has been only partially realized.
The interpretation of RS signals for forest canopies requires profound knowledge of the factors affecting their optical properties, which may be internal or external to the forest stand [6].To extract useful information related to canopy growth using time-series data, anomalies in time-series data that are unrelated to real changes in canopies should be eliminated, such as clouds and aerosol contamination [7,8] as well as bidirectional reflectance distribution function (BRDF) effects caused by topography and sun-sensor geometry variation [9][10][11][12][13][14]. Abundant studies have confirmed the existence of significant seasonal patterns in the optical RS signal trajectory that remain after removing the impacts of these factors external to the forest stand [11,13,15,16].Thus, variations in reflectance trajectories contain useful information related to factors internal to the stand.Based on previously reported experimental and modeling data, vegetation reflectance is primarily a function of tissue optical properties (reflectance and transmittance), canopy biophysical attributes (e.g., leaf area, foliage clumping) and background reflectance [17][18][19][20][21]. Leaf optical properties (LOPs) and leaf area index (LAI) are two of the main recognized internal factors involved in controlling canopy reflectance.
Seasonality in canopy spectral signals has been attributed to a varying LAI along with new leaf development and defoliation [22].In deciduous forests, LAI shows high seasonality and might be the most significant factor affecting canopy reflectance, as all leaves are shed in winter.However, the situation is very different for mature evergreen forests, which show a relatively stable total leaf area every year.Thus, changes in canopy reflectance do not necessarily imply changes in LAI [23].In evergreen forests, leaves have more than a one-year lifespan, and current-year new leaves remain throughout the winter.Although old leaves are shed every year, new leaves are also produced every year.The total leaf area in evergreen forests remains relatively stable throughout a year compared with deciduous forests.However, we still find significant seasonal variation in canopy reflectance in evergreen forests, which may not be caused by variations in LAI but rather by other internal factors.Further studies are needed to develop a more profound interpretation of the seasonality of optical RS signals for evergreen forests.
In addition to leaf area, LOPs are another significant factor that can change with time, which may strongly affect canopy reflectance [19,23,24], and leaf-age effects on canopy signals require more attention.Leaves affects the radiation field through its LOPs, including leaf reflectance and transmittance characteristics, which are wavelength-dependent [25,26].A few supporting studies have confirmed the combined effects of LAI and LOPs on the seasonality of gross ecosystem CO 2 exchange (GEE), photosynthesis and NIR reflectance for Amazonian forests [23,27,28].Many efforts aimed at accounting for effect of stand age on canopy reflectance with forest succession have improved the interpretation of canopy signals [29][30][31][32][33].However, only a few studies have indicated that in addition to stand age, leaf age might also have a significant impact.Our ability to both interpret RS signals and develop new RS technologies for vegetation depends directly on our ability to resolve the multitude of factors controlling canopy and landscape reflectance signatures.This situation inspired us to further evaluate robust biophysical interpretations of the seasonality of optical RS signals, with a focus on evergreen forests, by examining the effects of age-dependent leaf properties on canopy reflectance.
Different age cohorts of leaves coexist in the canopy of evergreen forests, which can be classified to two main age groups: current-year new leaves (≤1 a) and mature leaves (>1 a), which might exhibit different LOPs and, thus, different impacts on canopy reflectance.Mature leaves represent the majority and new leaves are the minority of leaves flushed every year.On the other hand, new leaves are mainly distributed at the top and in the outermost parts in tree crowns, while mature leaves are distributed in the lower and inner parts of the canopy.Thus, the seasonal variation of canopy reflectance might be strongly affected by changes in LOPs of both new leaves and mature leaves.This brings us to the crux of our study: the quantitative analysis and interpretation of the different leaf-aging effects of new leaves and mature leaves on seasonal variations in canopy reflectance.Hence, we addressed the following research questions in this study: (1) how do the LOPs of new leaves and mature leaves vary during the leaf maturation process?(2) How should the contributions of new leaves and mature leaves to canopy optical properties be quantified?(3) How does leaf aging affect the seasonal variation of the remotely sensed response spectral signals of evergreen forest canopies?Two main permanent sampling plots, ZH1 and WS3, were selected from S1 and S2, respectively, to measure and evaluate the effects of leaf aging on canopy optical properties.Two additional auxiliary plots, i.e., FZ1 and WS2, were also selected for comparison with the main plots.FZ1 exhibits the same plot conditions as ZH1 but is smaller, with a size of 30 × 40 m.WS3 is located next to WS2 and is eight years older than WS2.
Cunninghamia (Chinese fir) is a fast growing species, with its height increasing as much as 1 m per year [34], and usually matures approximately 20 a.The same Cunninghamia seedlings were planted in ZH1 and FZ1 after clear cutting.The stem density is the same in ZH1/FZ1, which exhibits half the density in WS2 (1920 trees ha −1 ) and WS3 (1967 trees ha −1 ).ZH1/FZ1 were planted with an initial stem density of 2500 trees ha −1 and thinned twice, once in 1997 and again in 2003, to a stable density of 1035 trees ha −1 .

Study Sites
This study focused on two National Research Stations of Forest Ecosystems in Huitong county, Hunan province, China, as shown in Figure 1.The first station (Station 1, S1) is located at 26°40′N, 109°26′E, and the second station (Station 2, S2) is located at 26°50′N, 109°45′E.S1 was established in 1983, and S2 was established in 1996.Two main permanent sampling plots, ZH1 and WS3, were selected from S1 and S2, respectively, to measure and evaluate the effects of leaf aging on canopy optical properties.Two additional auxiliary plots, i.e., FZ1 and WS2, were also selected for comparison with the main plots.FZ1 exhibits the same plot conditions as ZH1 but is smaller, with a size of 30 × 40 m.WS3 is located next to WS2 and is eight years older than WS2.
Cunninghamia (Chinese fir) is a fast growing species, with its height increasing as much as 1 m per year [34], and usually matures approximately 20 a.The same Cunninghamia seedlings were planted in ZH1 and FZ1 after clear cutting.The stem density is the same in ZH1/FZ1, which exhibits half the density in WS2 (1920 trees ha −1 ) and WS3 (1967 trees ha −1 ).ZH1/FZ1 were planted with an initial stem density of 2500 trees ha −1 and thinned twice, once in 1997 and again in 2003, to a stable density of 1035 trees ha −1 .
Figure 1.These four study plots are covered by Cunninghamia lanceolata (also known as Chinese fir) plantations, which were replanted after clear-cutting.

Canopy Structural Parameter Measurements
Crown shape measurements were taken from a total of forty trees in S2.The size of the selected trees was evenly distributed in terms of height (H) (from 5.4 m to 20 m) and diameter at breast height (DBH) (from 7.7 cm to 37.8 cm).Among the 40 trees, seven trees were located in plots close to the study site, and 33 trees were located within sites ZH1 and FZ1.The parameters measured for the characterization of crown shape included the following: crown width in the north-south direction

Canopy Structural Parameter Measurements
Crown shape measurements were taken from a total of forty trees in S2.The size of the selected trees was evenly distributed in terms of height (H) (from 5.4 m to 20 m) and diameter at breast height (DBH) (from 7.7 cm to 37.8 cm).Among the 40 trees, seven trees were located in plots close to the study site, and 33 trees were located within sites ZH1 and FZ1.The parameters measured for the characterization of crown shape included the following: crown width in the north-south direction (R1) and the east-west direction (R2), tree height (H1) and height under the crown (H2), from which we can obtain the height of the crown center (h).A detailed description of the crown shape measurements can be found in Appendix A.
Only DBH and H1 were measured annually for each tree in the study plots, and other canopy structure parameters were only measured once for the 40 trees.Hence, to characterize the dominant stand canopy structure changes with time, we needed to build allometric relationships between other unknown canopy structural parameters with DBH or H1.The derived canopy structure parameters included the following: Finally, the variations in canopy structure with stand development can be characterized by the means and standard deviations of the DBH and tree height (H1) for every individual tree.

LAI Measurements and Data Processing
Digital hemispherical photography (DHP) was the primary method for conducting regular monthly LAI measurements from 2005 until the present.From 2005 to 2006, photographs were taken by a worker every month using a CI-110 Plant Canopy Analyzer (Camas, WA, USA) to estimate the LAI.From 2007 onward; photographs were taken using a Canon EOS 40D digital camera equipped with a Nikon AF-DX 10.5 mm f/2.8 G ED fisheye lens.The camera was horizontally mounted at a fixed height of 0.2 m above the ground.The photographs were taken with automatic exposure under diffuse light conditions, typically soon before sunrise or after sunset.
In the ZH1 and FZ1 plots, measurements were taken on the 15th day of each month at five fixed locations per plot, facing in four cardinal directions.The images were processed using Gap Light Analyzer 2.0 software to calculate LAI DHP .LAI DHP measurements obtained with automatic exposure resulted in considerable underestimation because of underestimation of the ratio of green leaves to sky [35][36][37][38][39].We used the effective LAI (LAI e ) measured by two LAI-2000 Plant Canopy Analyzers (LAI-2000, LI-COR) to correct the system bias in LAI DHP data.One LAI-2000 unit was set at an open and flat area to measure diffuse sky light, and the other LAI-2000 unit was operated under the canopy at fixed locations as well as over the whole plots.
where ε can be considered as a systematic bias in the DHP method due to automatic exposure problems, and ε was set to 1.83 for the DHP method using the average value.The clumping index (Ω) was measured with the Tracing Radiation and Architecture of Canopies (TRAC, Natural Resources Canada) system to convert LAI e to true LAI (LAI t ) using the following equation: where the needle-to-shoot area ratio (γ e ) was set to 1.1 according to the results of the destructive sampling method described below, conducted in August 2015.The clumping index (Ω) was set to 0.8 and the value of the woody-to-total area ratio (α) was derived from the destructive samples for biomass estimation (0.2).
A detailed description of the data preprocessing methods of LAI can be found in Appendix B.

Spectral Measurements: Leaf and Soil
Soil and leaf samples were collected from WS2 and WS3 at site 2. Current-year leaves were too short to be measured in April.Therefore, we collected current-year shoots every month during the leaf expansion period from May to August 2017 to measure the seasonal variations in LOPs.Trees were selected from WS2 and WS3 separately to obtain branches from different age cohorts.Twenty branches were randomly selected and destructively harvested for each age group to conduct spectral measurements.Young leaves and mature leaves of 0-3 a were collected from the same branches; the leaves were carefully stored in sealed plastic bags and measured within 48 h.
Spectral measurements over the full spectral range (350-2500 nm) were carried out in the laboratory using a portable spectroradiometer (SVC HR-1024) attached to an integrating sphere (Model 1800-12S, Licor).All spectra were standardized using a barium sulfate standard and the calibrated light source supplied by Licor with the integrating sphere.We arranged 8-10 flat needles closely together for each measurement, avoiding overlaps or gaps, to obtain leaf samples that were sufficiently wide to cover the gap on integrating sphere.Wide transparent tape was used to assemble the leaf leaves on the vacuum side of the blade, which was removed from the central region to avoid any impacts on LOPs.LOPs were measured on the upper surface of the Cunninghamia lanceolata leaf samples.A group of leaf sample is shown in Figure A4 in Appendix A.

Landsat Observations
Landsat sensors have a long history and provide data with a fine spatial resolution (30 × 30 m) that can effectively capture forest stands in the landscape.In addition to their ideal spatial resolution, Landsat data provide certain other advantages, such as reducing the inconsistencies in observations by always viewing from almost the same direction (nadir view) and at the same time of day.However, one of the drawbacks of Landsat data (with a 16-day revisit period) is their relatively low temporal frequency, which is exacerbated by cloud-cover [40].We collected all available L1T Landsat TM/ETM+ images for our study site (path/row: 125/41) from 1987 to 2016, and abstracted pixel reflectance values with no cloud contamination.Given the high frequency of cloud cover in forested areas, Landsat observations over one year are insufficient to describe seasonal patterns.Therefore, all available pixel reflectance values from all years were sorted by the day of year (DOY), based on which seasonal trajectories of canopy reflectance were constructed.Landsat red and NIR band surface reflectance data were used to calculate EVI2, which takes advantage of the auto-correlative properties of surface reflectance spectra between the red and blue bands: where R NIR is reflectance in the near infrared band and R RED is reflectance in the red band.

Geometrical Optical Radiative Transfer (GORT) Model
The interaction between electromagnetic radiation and terrestrial plant canopies is a complex phenomenon and a key element in many RS applications.Among numerous methods for estimating the reflectance of forest canopies, the GORT model is based on the physical structure of the underlying scene.The GORT model is a hybrid of geometric optical (GO) and radiative transfer (RT) approaches for modeling canopy reflectance [41,42].The GO model [43][44][45] quantifies single scattering in the canopy well and captures the fundamental properties of the canopy bidirectional reflectance distribution function (BRDF).The assumptions of the GO model are that the scene is composed of three-dimensional solid objects on a contrasting background and that the overall canopy reflectance can be modeled as a weighted sum of the spectral signatures of its individual scene components, based on their corresponding areal proportions within a pixel.The RT model is used to describe the multiple scattering within canopy elements in the GORT model, and the GO model and RT model are linked using canopy gap probabilities [46,47].Due to the explicit consideration of crown gaps and mutual shadowing effects, the GORT model is suitable for simulating forest canopy reflectance with varying degrees of discontinuity.The GORT model has been successfully applied to predict the fundamental features of black spruce forest canopies [42], radiation penetrating the forest canopy to the forest floor [48] and spectral temporal manifestations of forest succession [49].
In the GO model [43][44][45], the reflectance of a pixel is modeled as the weighted sum of the spectral signatures of four scene components: sunlit ground, sunlit crown, shaded ground and shaded crown, as illustrated below: where S is the average reflectance of the forest canopy inside a pixel; and K g , K c , K z and K t represent the areal proportions of the four components in the pixel; thus, K g , K c , K z and K t sum to unity.G, C, Z and T are the corresponding spectral signatures (reflectance in a given wavelength range) of the four components, as shown in Figure 2, which are functions of the proportions of incoming directional beams and diffuse solar radiation.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 28 gaps and mutual shadowing effects, the GORT model is suitable for simulating forest canopy reflectance with varying degrees of discontinuity.The GORT model has been successfully applied to predict the fundamental features of black spruce forest canopies [42], radiation penetrating the forest canopy to the forest floor [48] and spectral temporal manifestations of forest succession [49].
In the GO model [43][44][45], the reflectance of a pixel is modeled as the weighted sum of the spectral signatures of four scene components: sunlit ground, sunlit crown, shaded ground and shaded crown, as illustrated below: where S is the average reflectance of the forest canopy inside a pixel; and , , and represent the areal proportions of the four components in the pixel; thus, , , and sum to unity., , and are the corresponding spectral signatures (reflectance in a given wavelength range) of the four components, as shown in Figure 2, which are functions of the proportions of incoming directional beams and diffuse solar radiation.The canopy structure parameters listed in Table 1 for the stands are derived based on the field measurements of the DBH and tree height of every individual tree in the stands, and the allometric relationships between DBH or tree height.Parameter h1 and h2 must represent the structure of the dominant canopy layer.Therefore, we used the mean height minus s.d. to derive h1 and the mean height plus s.d. to derive h2.FAVD is calculated by dividing LAI by the crown volume, by treating the tree crown as an ellipsoid.The canopy structure parameters listed in Table 1 for the stands are derived based on the field measurements of the DBH and tree height of every individual tree in the stands, and the allometric relationships between DBH or tree height.Parameter h1 and h2 must represent the structure of the dominant canopy layer.Therefore, we used the mean height minus s.d. to derive h1 and the mean height plus s.d. to derive h2.FAVD is calculated by dividing LAI by the crown volume, by treating the tree crown as an ellipsoid.

Contribution of Component LOPs to Canopy LOPs
Different leaf-age cohorts coexist in an individual tree canopy in evergreen forests and can be classified into two major groups: new leaves and mature leaves.The LOPs of new leaves are quite different from those of mature leaves, thus, neither of them is representative of the LOPs of the whole canopy.Therefore, we need to consider the contributions of new leaves and mature leaves together to obtain the overall canopy LOPs.

Leaf Area Proportion of New and Mature Leaves
The first factor affecting the average canopy spectral properties is the proportion of the leaf area at different ages: the higher the proportion of the leaf area for a given age, the greater its impact on canopy optical properties.Thus, we need to consider the proportions of new leaves and mature leaves viewed by the sensor.L total is the total LAI of all leaves, which consists of two parts: including the LAI of new leaves (L new ) and mature leaves (L mature ).The proportions of new and mature leaves were estimated based on destructive field measurements conducted by Zhongkun et al. [50], who studied the LAI of leaves at different leaf ages in Chinese fir, and we constructed the seasonally dynamic leaf area proportions by interpolating the data recorded using the phenology rule, assuming stable proportions to be reached at summer.

Spatial Organization of New and Mature Leaves
In addition to differences in the quantities of new and mature leaves, the spatial organization of new and mature leaves also has significant impacts on canopy component spectral signatures.Top crowns are mostly occupied by new leaves, and approximately 80% of the upper crown leaf area is occupied by new leaves in Chinese fir canopies [50].Leaves distributed in upper crowns have significant impacts on canopy reflectance when viewed from the top of the crown.In contrast, old leaves (1 a, 2 a, 3 a) are mainly located in the inner and lower parts of tree crowns, which are less likely to be observed directly from the nadir view.Thus, contribution of new leaves located at the top of canopies cannot be ignored.
When we simulate canopy reflectance based on Landsat viewing geometry, new leaves occupy the majority of the field of view of the sensor, and solar radiation interacts with the leaves in the top layer first before reaching the lower canopy.Thus, new leaves in the upper canopy have a larger influence on canopy reflectance than mature leaves in the lower canopy.The influences of new and mature leaves on canopy LOPs are modeled as the areal-weighted averages of new and mature leaves observed by the sensors, as follows: and where R ave and T ave are the average leaf reflectance and transmittance at the canopy scale, respectively.R new and R mature are the reflectance of new and mature leaves, respectively.T new and T mature are the transmittance of new and mature leaves, respectively.Finally, the parameters of w1 and w2 are the areal weights for new and mature leaves, respectively.
Here, we model the areal weights, w1 and w2, considering both leaf-area and leaf-age impacts, as follows: and where w1 highlights the importance of new leaves in the crown in the top canopy, and the weight of mature leaves (w2) is estimated from the total contribution of all leaves after subtracting the occlusion effect of the new leaves.

LOPs at the Canopy Scale
LOPs at the canopy scale are a combination LOPs of both new leaves and mature leaves.At the leaf scale, LOPs are age-dependent and the LOPs of new leaves and mature leaves vary differentially.At the canopy scale, age-dependent LOPs provide a mechanism for producing seasonally varying forest albedo and changing NIR to red ratios, independent of changes in other canopy attributes.Our hypothesis is that component LOPs vary with leaf maturation, and contribute to the seasonality in canopy reflectance trajectories.In the following part, we describe methods for retrieving seasonal LOPs and examining the relationship between LOPs and canopy signals.

Model Sensitivity Analysis
To retrieve LOPs from Landsat observations using the GORT model, we need to obtain a comprehensive understanding of the sensitivity of model driven parameters.First, wide ranges are allowed for all parameters, and we used canopy structural measurements of young stands (stand age = 1 a) as the lower limit and canopy structural parameters of mature stands (stand age = 33 a) as the upper limit.Stem density and background reflectance were set to their true values.Since sensitivity analysis covered parameter ranges on a long temporal scale (1 a to 33 a), insensitive parameters are eliminated.Then, the remaining sensitive parameters are involved in the new sensitivity analysis with a smaller range, taking the field measurements of each parameter as prior knowledge.
The global sensitivity of the model parameters is analyzed using the extension of the Fourier amplitude sensitivity testing (EFAST) method [51].EFAST is a variance decomposition method determining what fraction of the variance in the model output can be explained by the variation in each input parameter (i.e., partial variance).The basis of the EFAST method is a parametric transformation that can reduce multidimensional integrals over the input parametric space to one-dimensional quadratures using a search curve that scans the whole input space [52].Scanning is conducted so that each axis of the parametric space is explored at a different frequency.Then, Fourier decomposition is used to calculate both the first-order sensitivity (the contribution to the variance of the model output by each input, S i ) and total-order sensitivity (the first-order effect plus interactions with other inputs, S Ti ) of each input parameter, given as (Saltelli et al., 2008): and where X ∼i denotes the variation in all input parameters except for X i , and S ij is the contribution to the total variance from the interactions between parameters.Following Saltelli et al. [53], to compute S i and S Ti , we created a quasi-random sequence parameter sampling matrix, P, with dimensions of (m, n) for each SA test, where m is the sample size, and n is the number of input parameters.We set m = 2 n , which was sufficient to test the convergence of the sensitivity index and the stability of the rankings.Each row in matrix P represents a possible value set of X, and the quasi-random sequence helps to distribute the sampling points as uniformly as possible in the parameter space and avoid clustering, in addition to increasing the convergence rate.The global sensitivity analysis method is complex, and more details can be found in the work by Saltelli [53][54][55].Fortunately, SimLab software [40] can help implement the EFAST method.

Model Inversion Strategy
In this section, we describe the GORT model inversion strategy used to retrieve LOPs from satellite observations.Landsat reflectance in the red and NIR bands and the derived EVI2 were used to retrieve sensitive parameters.The multi-stage inversion strategy proposed by Li et al. [56] was adopted for parameter retrieval.The main objective of the iterative inversion process is to adjust the model parameters so that the model output reflectance is as close as possible to the observed reflectance.The most sensitive parameter is retrieved first, and used as prior knowledge in the next inversion stage.
One of the most popular methods for solving the inversion problem is to minimize the cost function of control variables.In this study, the cost function [57] used to retrieve sensitive parameters from Landsat surface reflectance data is as follows: where y n obs and f n (X) are the observed reflectance/EVI2 value and the corresponding modeled reflectance/EVI2 value, respectively.The variables δ n 2 and δ l 2 are the variances of the observational data and the prior distribution of parameters, respectively.The variables X l and X l prior are the parameter values and the initial values in the model, respectively.N is the number of observations, and L is the number of parameters.Sequential quadratic programming [58], an optimization algorithm for solving nonlinear programming problems, was adopted to search for the cost function minimum.

Validation: Direct and Indirect Methods
During the retrieval process, the LOPs of mature leaves in the red and NIR bands were set to fixed values, except for leaf transmittance in the NIR band, which was interpolated using field measurements; the LOPs of new leaves in the red band were also interpolated using SVC measurements for leaves.The retrieval results for new leaf LOPs in the NIR band were tested in direct and indirect ways, as shown in Figure 3.

Model Inversion Strategy
In this section, we describe the GORT model inversion strategy used to retrieve LOPs from satellite observations.Landsat reflectance in the red and NIR bands and the derived EVI2 were used to retrieve sensitive parameters.The multi-stage inversion strategy proposed by Li et al. [56] was adopted for parameter retrieval.The main objective of the iterative inversion process is to adjust the model parameters so that the model output reflectance is as close as possible to the observed reflectance.The most sensitive parameter is retrieved first, and used as prior knowledge in the next inversion stage.
One of the most popular methods for solving the inversion problem is to minimize the cost function of control variables.In this study, the cost function [57] used to retrieve sensitive parameters from Landsat surface reflectance data is as follows: where and ( ) are the observed reflectance/EVI2 value and the corresponding modeled reflectance/EVI2 value, respectively.The variables and are the variances of the observational data and the prior distribution of parameters, respectively.The variables and are the parameter values and the initial values in the model, respectively.N is the number of observations, and L is the number of parameters.Sequential quadratic programming [58], an optimization algorithm for solving nonlinear programming problems, was adopted to search for the cost function minimum.

Validation: Direct and Indirect Methods
During the retrieval process, the LOPs of mature leaves in the red and NIR bands were set to fixed values, except for leaf transmittance in the NIR band, which was interpolated using field measurements; the LOPs of new leaves in the red band were also interpolated using SVC measurements for leaves.The retrieval results for new leaf LOPs in the NIR band were tested in direct and indirect ways, as shown in Figure 3.We decomposed retrieved LOPs into new leaf (time variant) and mature leaf (known) components according to the calculated the contribution weights (w1 and w2) described in Section 3.2, and evaluated the decomposed LOPs transmittance using SVC field measurements directly.Additionally, we tested the applicability of the retrieved LOPs as a second-step validation.LOPs retrieved at one site (ZH1 plot) were applied to another site (FZ1 plot) to simulate canopy reflectance using the GORT forward simulating mode, and evaluated pixel reflectance using Landsat observations.FZ1 plot and ZH1 plot have similar growing conditions, but are located at different stands at Station 1.These sites exhibit different canopy structure properties but with the same tree species, stem density and stand age and close LAI values.In the simulation of canopy reflectance for FZ1 plot, the same LOPs, sun-sensor geometry and soil reflectance for the ZH1 plot were used, but a used different canopy structure was employed.
During the forward simulation process for canopy reflectance, three situations were considered and compared: (1) without taking leaf-age effects into consideration: setting LOPs parameters as fixed values using monthly average r L -nir and t L -nir of mature leaves; (2) considering the leaf-age effects for mature leaves: setting LOPs parameters using monthly varying r L -nir and t L -nir datasets of mature leaves; and (3) considering the leaf-age effects of both mature leaves and new leaves: setting LOPs parameters using monthly varying canopy-scale r L -nir and t L -nir datasets, which are up-scaled from the seasonal LOPs of new leaves and mature leaves.Finally, simulated canopy reflectance signatures were compared with pixel reflectance derived from Landsat time-series observations at FZ1 as an indirect validation.

Total-Order and Single-Order Sensitivity Analysis Results
Sensitivity analyses run two times for the GORT model for the red band and the NIR band.In the first sensitivity analysis, parameters with wide ranging variations were involved in the sensitivity analysis.In the second sensitivity analysis, we eliminated stem density (λ) and crown radius (r) in the analysis and focused on the remaining parameters, since stem density was known during our study time, and the average tree crown radius can be estimated based on the allometric relationship with DBH.In the second sensitivity analysis, the remaining parameters (LAI, h1, h2, λ, r L , t L , r G ) varied within the same range as in the first sensitivity analysis.
The sensitivity of the model input parameters was characterized using the total-order sensitivity index (Figure 4).For canopy reflectance in the red band, canopy reflectance is most sensitive to the crown radius and stem density in the first step of sensitivity analysis (Figure 4A).Given r and λ are known, LAI, h2, r L and t L become the most influential parameters in order (Figure 4B).For canopy reflectance in the NIR band, r L , h2 and r are the main influential parameters over the long-term (Figure 4A).After removing the uncertainty of the crown radius and stem density, t L becomes the third most influential parameter following r L and h2.
To reveal the potential for retrieving uncertain parameters from canopy reflectance, we further analyzed the single-order sensitivity of each sensitive parameter in the GORT model outputs by fixing other parameters at stand average values.Taking the year of 2005 as an example, the single-order sensitivity results are shown in Figure 5.We found that during our study period, LAI and h2 were not influential regarding canopy spectral signatures in the red (Figure 5A) band, since LAI remains within a limited range (LAI > 3 m 2 /m 2 ) for mature evergreen forests, and the crown radius remained the most influential parameter, especially in the red band (Figure 5A).For canopy reflectance in the NIR band, h2 was not influential within the range (18-23 m) for mature evergreen forests (Figure 5B).LOPs were sensitive parameters for both the red and NIR bands.

Prior Knowledge of Model Parameters
All prior knowledge of parameter values was obtained from field measurements and used for analyzing the sensitivity indexes of model parameters, as summarized in Table 2.According to totalorder and first-order sensitivity analysis results, the influential parameters include: the crown radius (r) , and , which were assumed to be adjustable during retrieval process.Non-influential parameters were set to values estimated using field measurements as described in Section 2.2.Multiprocess model inversion was conducted in the following sequences: (1) Landsat red band reflectance −> crown radius (r); (2) Landsat NIR band reflectance/EVI2 −> leaf reflectance ( ) for the NIR band; (3) Landsat NIR band reflectance/EVI2 −> leaf transmittance ( ) for the NIR band; (4) Landsat red band reflectance/EVI2 −> leaf reflectance ( ) for the red band; (5) Landsat red band reflectance/EVI2 −> leaf transmittance ( ) for the red band.The results for the retrieved parameters are summarized in Table 2.

Prior Knowledge of Model Parameters
All prior knowledge of parameter values was obtained from field measurements and used for analyzing the sensitivity indexes of model parameters, as summarized in Table 2.According to totalorder and first-order sensitivity analysis results, the influential parameters include: the crown radius (r) , and , which were assumed to be adjustable during retrieval process.Non-influential parameters were set to values estimated using field measurements as described in Section 2.2.Multiprocess model inversion was conducted in the following sequences: (1) Landsat red band reflectance −> crown radius (r); (2) Landsat NIR band reflectance/EVI2 −> leaf reflectance ( ) for the NIR band; (3) Landsat NIR band reflectance/EVI2 −> leaf transmittance ( ) for the NIR band; (4) Landsat red band reflectance/EVI2 −> leaf reflectance ( ) for the red band; (5) Landsat red band reflectance/EVI2 −> leaf transmittance ( ) for the red band.The results for the retrieved parameters are summarized in Table 2.

Prior Knowledge of Model Parameters
All prior knowledge of parameter values was obtained from field measurements and used for analyzing the sensitivity indexes of model parameters, as summarized in Table 2.According to total-order and first-order sensitivity analysis results, the influential parameters include: the crown radius (r), r L and t L , which were assumed to be adjustable during retrieval process.Non-influential parameters were set to values estimated using field measurements as described in Section 2.2.Multi-process model inversion was conducted in the following sequences: (1) Landsat red band reflectance −> crown radius (r); (2) Landsat NIR band reflectance/EVI2 −> leaf reflectance (r L ) for the NIR band; (3) Landsat NIR band reflectance/EVI2 −> leaf transmittance (t L ) for the NIR band; (4) Landsat red band reflectance/EVI2 −> leaf reflectance (r L ) for the red band; (5) Landsat red band reflectance/EVI2 −> leaf transmittance (t L ) for the red band.The results for the retrieved parameters are summarized in Table 2.We first studied the leaf-age effect on leaf optical properties, and found that leaves in the canopy can be classified into two age groups: new leaves (0 a) and mature leaves (1-3 a).Field SVC measurements of the full spectra for all leaf samples were converted to Landsat-view using the Landsat relative spectral response (RSR) functions in the following three bands: green, red and NIR.Leaf samples were grouped into four age classes: 0 a, 1 a, 2 a and 3 a.Significant differences in LOPs were observed between 0 and 1 a leaves, while the differences between mature leaves at different ages (1-3 a) were not distinct (Figure 6).Thus, LOPs mainly varied within 0-1 a, i.e., during the maturation process of newly flushed leaves.We first studied the leaf-age effect on leaf optical properties, and found that leaves in the canopy can be classified into two age groups: new leaves (0 a) and mature leaves (1-3 a).Field SVC measurements of the full spectra for all leaf samples were converted to Landsat-view using the Landsat relative spectral response (RSR) functions in the following three bands: green, red and NIR.Leaf samples were grouped into four age classes: 0 a, 1 a, 2 a and 3 a.Significant differences in LOPs were observed between 0 and 1 a leaves, while the differences between mature leaves at different ages (1-3 a) were not distinct (Figure 6).Thus, LOPs mainly varied within 0-1 a, i.e., during the maturation process of newly flushed leaves.Differences in leaf reflectance at different ages are observed in two visible bands and the NIR band (Figure 6A1-C1).New Leaves (0 a) exhibited a higher reflectance in the green and red bands compared with mature leaves (Figure 6A1-B1).In the NIR band, 1 a and 2 a leaves exhibited a slight increase in reflectance, but the reflectance of 3 a leaves was similar to that of 0 a leaves (Figure 6C1).
As can be observed in the Landsat-view transmittance signatures shown in Figure 6A2-C2, the trends in the changes in transmittance characteristics as a function of age were similar for these three bands.Significant differences were observed between 0 a and 1 a leaves, with the 0 a leaves showing Differences in leaf reflectance at different ages are observed in two visible bands and the NIR band (Figure 6A1-C1).New Leaves (0 a) exhibited a higher reflectance in the green and red bands compared with mature leaves (Figure 6A1-B1).In the NIR band, 1 a and 2 a leaves exhibited a slight increase in reflectance, but the reflectance of 3 a leaves was similar to that of 0 a leaves (Figure 6C1).
As can be observed in the Landsat-view transmittance signatures shown in Figure 6A2-C2, the trends in the changes in transmittance characteristics as a function of age were similar for these three bands.Significant differences were observed between 0 a and 1 a leaves, with the 0 a leaves showing a consistently higher level of transmittance compared with mature leaves.The decrease in transmittance as a function of increasing leaf age is due to the increase in absorption characteristics after a new leaf matures, which is also supported by the increasing area between the upper and lower set of curves shown in Figure 7A1-E1.
Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 28 a consistently higher level of transmittance compared with mature leaves.The decrease in transmittance as a function of increasing leaf age is due to the increase in absorption characteristics after a new leaf matures, which is also supported by the increasing area between the upper and lower set of curves shown in Figure 7A1-E1.

Leaves in Different Seasons
The results of monthly LOPs measurements during the leaf expansion period (May to October) of 0 a leaves sampled from the Chinese fir trees are presented in Figure 7.As shown in the SVC spectral curve data for Chinese fir, significant differences were observed in the following spectral regions: the green peak (≈530-600 nm), the chlorophyll absorption well (≈660-690 nm) and along the NIR plateau (750-900 nm).
Differences in leaf reflectance are most pronounced in the green band.New leaves (0 a) exhibited a consistently lower green peak reflectance from May to October, while reflectance in the chlorophyll absorption region was shown to increase with age, especially during the first two months of the leaf expansion period, which is consistent with the increasing trend in NIR reflectance from 0 a to 1 a leaves as presented in Figure 6C1.Although this variation in NIR band reflectance is not as pronounced as in the green bands, it could have a greater impact at the canopy scale from the view of satellites, as illustrated in the following section.
With regard to leaf-level transmittance (Figure 7A1-E1), 0 a leaves showed a remarkable decrease in transmittance in the visible region from May to October, which was similar to that observed in the variations in leaf transmittance from 0 a to 1 a (Figure 7A2-B2).NIR transmittance characteristics first showed a slight increasing trend (from 1 May to 28 July) and presented a decreasing trend thereafter (from 28 July to 13 October), which coincided with the changing trend in new leaf NIR transmittance trajectory retrieved from Landsat observations (Figure 8A1).However, the observations from May to October could not fully explain the gap in the NIR transmittance of 0 a and 1 a leaves (Figure 8C2).One possible reason is that leaf transmittance decreases during winter (November to April), which may be supported by the new leaf transmittance retrieved from Landsat observations in winter time.As can be observed in Figure 8C1, NIR transmittance decreased significantly in November and December, but there is still a gap between the retrieved new leaf Field spectral curve data (mean ± 1 s.d.) in shortwave bands for new (A1-E1) and mature (A2-E2) needle samples (n) from Chinese fir collected on 5 May, 24 June, 28 July, 14 September and 13 October 2017.From May to October, n = 22, 15, 26, 30 and 39 respectively, and each sample includes 5 to 8 leaves.Reflectance data (r L ) for each month are presented as the lower set of curves within each plot, while transmittance data (t L ) are presented as the upper set of curves.Absorption characteristics are depicted based on the area between the upper and lower set of curves.SD values for new leaf t L are depicted in the upper and lower dash lines, to avoid overlap with r L , while SD values for all other r L and t L measurements are depicted in gray buffed areas.

Leaves in Different Seasons
The results of monthly LOPs measurements during the leaf expansion period (May to October) of 0 a leaves sampled from the Chinese fir trees are presented in Figure 7.As shown in the SVC spectral curve data for Chinese fir, significant differences were observed in the following spectral regions: the green peak (≈530-600 nm), the chlorophyll absorption well (≈660-690 nm) and along the NIR plateau (750-900 nm).
Differences in leaf reflectance are most pronounced in the green band.New leaves (0 a) exhibited a consistently lower green peak reflectance from May to October, while reflectance in the chlorophyll absorption region was shown to increase with age, especially during the first two months of the leaf expansion period, which is consistent with the increasing trend in NIR reflectance from 0 a to 1 a leaves as presented in Figure 6C1.Although this variation in NIR band reflectance is not as pronounced as in the green bands, it could have a greater impact at the canopy scale from the view of satellites, as illustrated in the following section.
With regard to leaf-level transmittance (Figure 7A1-E1), 0 a leaves showed a remarkable decrease in transmittance in the visible region from May to October, which was similar to that observed in the variations in leaf transmittance from 0 a to 1 a (Figure 7A2-B2).NIR transmittance characteristics first showed a slight increasing trend (from 1 May to 28 July) and presented a decreasing trend thereafter (from 28 July to 13 October), which coincided with the changing trend in new leaf NIR transmittance trajectory retrieved from Landsat observations (Figure 8A1).However, the observations from May to October could not fully explain the gap in the NIR transmittance of 0 a and 1 a leaves (Figure 8C2).One possible reason is that leaf transmittance decreases during winter (November to April), which may be supported by the new leaf transmittance retrieved from Landsat observations in winter time.As can be observed in Figure 8C1, NIR transmittance decreased significantly in November and December, but there is still a gap between the retrieved new leaf transmittance (0.2) and measured leaf transmittance (0.25) at the point of reaching 1 a in May of the following year after leaf production.
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 28 transmittance (0.2) and measured leaf transmittance (0.25) at the point of reaching 1 a in May of the following year after leaf production.

Leaf-Age Effects on Variability in Landsat-Viewed Canopy Reflectance
Landsat records canopy LOPs from the sensor view and time-series data can be used to estimate seasonal variations in LOPs.First, we attempted to retrieve LOPs at the canopy scale from Landsat observations using the GORT model.Then, we estimated the new leaf component from the retrieved results by setting mature leaf component parameters as known parameters using field measurements.The results were evaluated via direct and indirect methods.

Leaf Optical Properties at the Canopy Scale
After considering the uncertainties and sensitivities of other parameters in the GORT model, seasonal LOPs were retrieved from Landsat observations using the proposed multi-stage inversion method by minimizing the cost function.Retrieved LOPs at the canopy scale are shown in Figure 8 and compared with field measurements of the spectral signatures of leaves.Changes in Landsat reflectance were associated with variations in leaf optical properties.As shown in Figure 8A1-B1, the most dramatic changes occurred in the NIR band, resulting in increased reflectance and transmittance during spring (April to September) and decreased reflectance and transmittance during winter (October to December).However, variations in NIR reflectance (0.5 to 0.55) were less notable than that of NIR transmittance (0.25 to 0.4).In the red band (Figure 8A2-B2), differences in LOPs mainly occurred in the first three months (May to July) and then became stable in the rest of the year.Both red reflectance and transmittance presented low values, and transmittance was slightly lower than reflectance.No LOPs changes in red band could be observed from Landsat alone, but including EVI2 data helped to some extent.

Seasonal Leaf Optical Properties
When viewed from the top of the crown, LOPs retrieved from Landsat observations are a combination of leaves of all ages.As analyzed in Section 3.3, we account for both the effects of leaf area and its spatial organization in the canopy to estimate the contribution of different leaves.Figure 9 shows the average seasonal variations in total leaf area (A) and leaf proportions (B) of new leaves and mature leaves in the canopy, and their contributions to canopy spectral properties (C).The time

Leaf-Age Effects on Variability in Landsat-Viewed Canopy Reflectance
Landsat records canopy LOPs from the sensor view and time-series data can be used to estimate seasonal variations in LOPs.First, we attempted to retrieve LOPs at the canopy scale from Landsat observations using the GORT model.Then, we estimated the new leaf component from the retrieved results by setting mature leaf component parameters as known parameters using field measurements.The results were evaluated via direct and indirect methods.

Leaf Optical Properties at the Canopy Scale
After considering the uncertainties and sensitivities of other parameters in the GORT model, seasonal LOPs were retrieved from Landsat observations using the proposed multi-stage inversion method by minimizing the cost function.Retrieved LOPs at the canopy scale are shown in Figure 8 and compared with field measurements of the spectral signatures of leaves.Changes in Landsat reflectance were associated with variations in leaf optical properties.As shown in Figure 8A1-B1, the most dramatic changes occurred in the NIR band, resulting in increased reflectance and transmittance during spring (April to September) and decreased reflectance and transmittance during winter (October to December).However, variations in NIR reflectance (0.5 to 0.55) were less notable than that of NIR transmittance (0.25 to 0.4).In the red band (Figure 8A2-B2), differences in LOPs mainly occurred in the first three months (May to July) and then became stable in the rest of the year.Both red reflectance and transmittance presented low values, and transmittance was slightly lower than reflectance.No LOPs changes in red band could be observed from Landsat alone, but including EVI2 data helped to some extent.

Seasonal Leaf Optical Properties
When viewed from the top of the crown, LOPs retrieved from Landsat observations are a combination of leaves of all ages.As analyzed in Section 3.3, we account for both the effects of leaf area and its spatial organization in the canopy to estimate the contribution of different leaves.Figure 9 shows the average seasonal variations in total leaf area (A) and leaf proportions (B) of new leaves and mature leaves in the canopy, and their contributions to canopy spectral properties (C).The time period starts at leaf expansion at age 0 (April) and ends with leaf maturity (March in the following year) at age 1, following the same twelve-month cycle with new leaf expansion each year.We can see that LAI does not change very dramatically for evergreen forests between seasons in our study period (Figure 9A).Thus, the satellite-observed seasonality of canopy reflectance is not determined by LAI, but rather is the result of variations in the LOPs of new leaves during the maturation process.Although new leaves only account for approximately 30% of the total leaf area (Figure 9B), new leaves made high contribution (approximately 80% at the peak time in summer) to canopy spectral properties from the crown top view (Figure 9C).
Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 28 period starts at leaf expansion at age 0 (April) and ends with leaf maturity (March in the following year) at age 1, following the same twelve-month cycle with new leaf expansion each year.We can see that LAI does not change very dramatically for evergreen forests between seasons in our study period (Figure 9A).Thus, the satellite-observed seasonality of canopy reflectance is not determined by LAI, but rather is the result of variations in the LOPs of new leaves during the maturation process.
Although new leaves only account for approximately 30% of the total leaf area (Figure 9B), new leaves made high contribution (approximately 80% at the peak time in summer) to canopy spectral properties from the crown top view (Figure 9C).The seasonal trajectories of mature leaf LOPs (Figure 10 B1,B2) were interpolated from field measurements, which were quite stable during growing season, with only one exception: NIR transmittance.Mature leaf transmittance in the NIR increased by 30% over the three months following the beginning of a new growing cycle in May (Figure 10B1), while mature leaves showed no significant difference in NIR reflectance during the whole growing season (May to November).
New leaf reflectance or transmittance can be estimated by deducting the contribution of mature leaves from the retrieved LOPs results, which are shown in Figure 10.It can be seen that the retrieved new leaf LOPs were in good agreement with the field observations (Figure 10A1,A2) from May to October.We can also see that new leaf LOPs in the red band decreased in the first three months (from April to June) and became stable thereafter (Figure 10A2).More importantly, new leaf NIR REF continued to increase from May to September after production (approximately 11%), while new leaf NIR TRA follow a nonlinear trajectory, increasing from May to August (approximately 16%) and decreasing thereafter (approximately 14%).A small increase in new leaf NIR REF and a small decrease in new leaf NIR TRA could both translate into a larger impact at the canopy scale, when the canopy is coated with new leaves, and the impact of the leaf-age effect at the canopy scale will be further studied in the following section.
A gap exists between the estimated new leaf NIR transmittance and SVC measurements in May, and the s.d.values of new leaf field measurements in May are greater than the measurements in other months (Figure 10A1,A2).This might be explained by the overestimation of measured leaf transmittance samples when leaves were small, which was caused by unavoidable small gaps because of minor misalignment in the arrangement of leaves.Large s.d.values for new leaf field measurements in May could also occur due to the differences in leaf expansion rate and maturity times in the early leaf flush period.The seasonal trajectories of mature leaf LOPs (Figure 10 B1,B2) were interpolated from field measurements, which were quite stable during growing season, with only one exception: NIR transmittance.Mature leaf transmittance in the NIR increased by 30% over the three months following the beginning of a new growing cycle in May (Figure 10B1), while mature leaves showed no significant difference in NIR reflectance during the whole growing season (May to November).
New leaf reflectance or transmittance can be estimated by deducting the contribution of mature leaves from the retrieved LOPs results, which are shown in Figure 10.It can be seen that the retrieved new leaf LOPs were in good agreement with the field observations (Figure 10A1,A2) from May to October.We can also see that new leaf LOPs in the red band decreased in the first three months (from April to June) and became stable thereafter (Figure 10A2).More importantly, new leaf NIR REF continued to increase from May to September after production (approximately 11%), while new leaf NIR TRA follow a nonlinear trajectory, increasing from May to August (approximately 16%) and decreasing thereafter (approximately 14%).A small increase in new leaf NIR REF and a small decrease in new leaf NIR TRA could both translate into a larger impact at the canopy scale, when the canopy is coated with new leaves, and the impact of the leaf-age effect at the canopy scale will be further studied in the following section.
A gap exists between the estimated new leaf NIR transmittance and SVC measurements in May, and the s.d.values of new leaf field measurements in May are greater than the measurements in other months (Figure 10A1,A2).This might be explained by the overestimation of measured leaf transmittance samples when leaves were small, which was caused by unavoidable small gaps because of minor misalignment in the arrangement of leaves.Large s.d.values for new leaf field measurements in May could also occur due to the differences in leaf expansion rate and maturity times in the early leaf flush period.

Leaf-Age Effects at Pixel Scale Based on Satellite Observations
LAI and sun-sensor geometry are widely recognized as the main factors contributing to seasonality in RS signals [12,19].Thus, we need to distinguish leaf-age effects from the contributions of these two factors.To identify the contribution of leaf-age effects to the seasonality in Landsat signals, we compared the differences in the results with/without considering leaf-age effects.The forward-mode GORT model was used to simulate canopy reflectance at the pixel scale with Landsatviewing geometry, and the model-driven parameters are listed in Table 1.We simulate canopy reflectance at NIR band and red band in the following three circumstances: 1.In the first circumstance, we ignore the leaf-age effects caused by aging mature leaves and growing new leaves, and only consider variations in LAI and sun geometry using field data.
LOPs in the GORT model were fixed using the mean LOPs for mature leaves measured from May to September in 2017 (NIR band: REFmature = 0.51,TRAmature = 0.26 or 0.34; Red band: REFmature = 0.06,TRAmature = 0.01).2. Based on circumstance 1, we include variations in LOPs caused by aging mature leaves.For comparative purposes, we first added variations in mature leaves using data shown in Figure 10B1 to drive the GORT model.3. Based on circumstance 2, we further include variations in LOPs caused by production and expansion of new leaves.LOPs for the GORT model are shown in Figure 8B1.LOPs at the canopy scale retrieved at the ZH1 site are applied to the FZ1 site with different canopy structure parameters.
The overall GORT-simulated results obtained for the three circumstances are compared in Figure 11.After considering the mature leaf-age effects on canopy reflectance, the R 2 between the simulated canopy reflectance/EVI2 and the Landsat observations increased significantly (from 0.38 to 0.83 for canopy NIR reflectance and from 0.26 to 0.43 for canopy EVI2), as shown in Figure 11.However, the R 2 of both groups in the GORT simulation results was poor in the red band and RMSE remained low.
In circumstance 1, we can see that LAI is stable (Figure 9A) in these sub-tropical forests and is not the main factor contributing to the seasonality of canopy NIR REF.Seasonal variation in sunsensor geometry from January to June to December (SZN varies from 56° to 25° to 56°) causes a small amount of seasonality of forest albedo, as shown based on the blue line in Figure 11A1.However, SZN alone does not sufficiently explain the seasonality of the Landsat canopy signals in our study

Leaf-Age Effects at Pixel Scale Based on Satellite Observations
LAI and sun-sensor geometry are widely recognized as the main factors contributing to seasonality in RS signals [12,19].Thus, we need to distinguish leaf-age effects from the contributions of these two factors.To identify the contribution of leaf-age effects to the seasonality in Landsat signals, we compared the differences in the results with/without considering leaf-age effects.The forward-mode GORT model was used to simulate canopy reflectance at the pixel scale with Landsat-viewing geometry, and the model-driven parameters are listed in Table 1.We simulate canopy reflectance at NIR band and red band in the following three circumstances: 1.
In the first circumstance, we ignore the leaf-age effects caused by aging mature leaves and growing new leaves, and only consider variations in LAI and sun geometry using field data.LOPs in the GORT model were fixed using the mean LOPs for mature leaves measured from May to September in 2017 (NIR band: REF mature = 0.51, TRA mature = 0.26 or 0.34; Red band: REF mature = 0.06, TRA mature = 0.01).

2.
Based on circumstance 1, we include variations in LOPs caused by aging mature leaves.
For comparative purposes, we first added variations in mature leaves using data shown in Figure 10B1 to drive the GORT model.

3.
Based on circumstance 2, we further include variations in LOPs caused by production and expansion of new leaves.LOPs for the GORT model are shown in Figure 8B1.LOPs at the canopy scale retrieved at the ZH1 site are applied to the FZ1 site with different canopy structure parameters.
The overall GORT-simulated results obtained for the three circumstances are compared in Figure 11.After considering the mature leaf-age effects on canopy reflectance, the R 2 between the simulated canopy reflectance/EVI2 and the Landsat observations increased significantly (from 0.38 to 0.83 for canopy NIR reflectance and from 0.26 to 0.43 for canopy EVI2), as shown in Figure 11.However, the R 2 of both groups in the GORT simulation results was poor in the red band and RMSE remained low.
In circumstance 1, we can see that LAI is stable (Figure 9A) in these sub-tropical forests and is not the main factor contributing to the seasonality of canopy NIR REF.Seasonal variation in sun-sensor geometry from January to June to December (SZN varies from 56 • to 25 • to 56 • ) causes a small amount of seasonality of forest albedo, as shown based on the blue line in Figure 11A1.However, SZN alone does not sufficiently explain the seasonality of the Landsat canopy signals in our study sites.Assuming that the LOPs of mature leaves do not change, if we apply mature leaf LOPs measured at the peak time in summer (NIR TRA = 0.34) to other times, we might overestimate canopy NIR REF in winter; and if we apply mature leaf LOPs measured in early spring (NIR TRA = 0.26), we might underestimate canopy NIR REF in summer.
In circumstance 2, we consider the leaf-aging effects of the mature leaves, which constitute the majority (>60%) of the canopy LAI (Figure 9B).The LOPs of mature leaves are relatively stable in the red and NIR bands, therefore, mature leaves have no impact on canopy red REF.However, there is one notable exception: mature leaf NIR TRA is lower in winter and greater in summer and continued to increase from May to October (Figure 7A2-E2).The winter-summer differences in leaf NIR TRA contribute to a steeper NIR REF/EVI2 trajectory from April to July at the canopy scale, as shown based on the area between the blue and gray curves in Figure 11A1,C1.
In circumstance 3, we can observe how variations in the LOPs of new leaves contributed to the seasonal variation in canopy NIR reflectance and derived EVI2 (Figure 11A1,C1).As shown based on the area between the green and gray curves in Figure 11A1,C1, EVI2 trajectories follow the same pattern (Figure 11C1) as NIR reflectance.In general, new leaf production and expansion increases canopy NIR reflectance from May to October (Figure 11A1).An increased new leaf TRA (18%) also contributes to the increase in NIR REF at the canopy scale, but only before August.After September, an increase in NIR REF is caused by an increase in leaf REF.However, the changing trend of canopy NIR REF is nonlinear, which increases from May to August and decreases from August to October.As illustrated in Figure 10A1, canopy NIR REF increases from May to August partially due to increasing new leaf NIR REF during this period; conversely, canopy NIR REF decreases from August to October due to decreasing new leaf TRA at this interval.
The upper-layer-located new leaves provide a mechanism for producing greater seasonality of forest albedo in addition to mature leaves.A small increase in new leaf NIR REF (0.05 unit) has significant impacts at the canopy scale, since new leaves are generally located at the top of the canopy.New leaf expansion also contributes to negligible increase in canopy Red REF from March to July (Figure 11B1).EVI2 trajectories follow the same patterns (Figure 11C1) as NIR reflectance trajectories, arising from the linear dependence of EVI2 on NIR reflectance, as proved by a previous study [23].We can conclude that both new and mature leaves contribute to the seasonality of forest albedo, which is independent of changes in other canopy attributes.In circumstance 2, we consider the leaf-aging effects of the mature leaves, which constitute the majority (>60%) of the canopy LAI (Figure 9B).The LOPs of mature leaves are relatively stable in the red and NIR bands, therefore, mature leaves have no impact on canopy red REF.However, there is one notable exception: mature leaf NIR TRA is lower in winter and greater in summer and continued to increase from May to October (Figure 7A2-E2).The winter-summer differences in leaf NIR TRA contribute to a steeper NIR REF/EVI2 trajectory from April to July at the canopy scale, as shown based on the area between the blue and gray curves in Figure 11A1,C1.
In circumstance 3, we can observe how variations in the LOPs of new leaves contributed to the seasonal variation in canopy NIR reflectance and derived EVI2 (Figure 11A1,C1).As shown based on the area between the green and gray curves in Figure 11A1,C1, EVI2 trajectories follow the same pattern (Figure 11C1) as NIR reflectance.In general, new leaf production and expansion increases canopy NIR reflectance from May to October (Figure 11A1).An increased new leaf TRA (18%) also contributes to the increase in NIR REF at the canopy scale, but only before August.After September, an increase in NIR REF is caused by an increase in leaf REF.However, the changing trend of canopy NIR REF is nonlinear, which increases from May to August and decreases from August to October.As illustrated in Figure 10A1, canopy NIR REF increases from May to August partially due to increasing new leaf NIR REF during this period; conversely, canopy NIR REF decreases from August to October due to decreasing new leaf TRA at this interval.
The upper-layer-located new leaves provide a mechanism for producing greater seasonality of forest albedo in addition to mature leaves.A small increase in new leaf NIR REF (0.05 unit) has significant impacts at the canopy scale, since new leaves are generally located at the top of the canopy.New leaf expansion also contributes to negligible increase in canopy Red REF from March to July (Figure 11B1).EVI2 trajectories follow the same patterns (Figure 11C1) as NIR reflectance trajectories, arising from the linear dependence of EVI2 on NIR reflectance, as proved by a previous study [23].We can conclude that both new and mature leaves contribute to the seasonality of forest albedo, which is independent of changes in other canopy attributes.

Spectral Changes and Leaf Aging
For Chinese fir, the leaf spectra mainly vary during the first year of leaf production (0-1 a), and then remain relatively stable from 1-3 a.During the new leaf expansion process, leaf water content (LWC) and the specific leaf area (SLA) decrease rapidly, resulting in variations in LOPs: including decreased transmittance and increased absorbance in the visible and NIR bands and decreased visible, but increased NIR reflectance.Leaf optical properties show strong age dependence and great seasonal difference for young and mature leaves.Spectral changes in the visible portion of the spectrum characterized the new leaves while increased NIR transmittance and decreased NIR absorbance characterized mature leaves.Changes observed in the visible spectrum matched similar observations for Carica [59] and Aldina [60].In general, absorbance at the green peak (550 nm) increases, while reflectance and transmittance decrease at this band, which could be attributed to an increase in chlorophyll [59] and changes in leaf internal anatomy [61].However, the changes of LOPs in the NIR band are complicated.
In agreement with previous findings [60], mature leaves were found to show a lower leaf albedo (leaf reflectance plus transmittance) at the NIR band than new leaves due to leaf aging.Early changes in the NIR band were noted during new leaf expansion, with a slight increase (11%) in leaf NIR REF being observed from May to October and a continuous increase (14%) in leaf TRA being recorded from May to August, followed by a slight decrease in NIR TRA during autumn and winter.Some studies [62] have found a dramatic increase in leaf NIR reflectance between July and August, which might be caused by differences in the leaf characteristics.Low-wax leaves display a continuous increase in NIR reflectance until maturity, while high-wax leaves show little increase in NIR reflectance [63].Epicuticular waxes and thick cuticles could mask the effect of NIR increases caused by inter-cellular development in young high-wax leaves.These characteristics might explain why only a slight increase in NIR reflectance was found during leaf maturation in the present study, since Chinese fir is a high-wax species.
In contrast to previous studies, we observed a continuous increase in NIR transmittance for mature leaves (1-3 a).However, a significant gap in NIR transmittance existed between new leaves (measured from May to October) and mature leave (measured in May), and we only observed a small decrease in NIR transmittance for new leaves from September to October.However, a dramatic decrease in NIR transmittance might happen in winter time (November to December to January) according to the results retrieved from Landsat observations, which will be further studied in the future.Similar observations have been made by others [60], who found a consistent decrease in NIR transmittance during the last nine months of the leaf cycle.

Leaf-Aging Effects on Canopy Reflectance
At the canopy scale, canopy reflectance is a product of competing mechanisms of light absorption [60].The interaction of photons with dense forests is characterized by strong scattering in the NIR and equally strong absorption in the shorter red and blue bands.The NIR reflectance of these forests is an order of magnitude greater than the reflectance at red (blue) band.On the other hand, aerosol scattering has greater impacts on reflectance in visible bands than the NIR band.In the NIR band, any mechanism that increases leaf absorption (such as decreased transmittance or reflectance) will have an enhanced effect on canopy reflectance [60].Thus, we could expect leaf aging to have its largest impact in the NIR band.
Age-related LOPs are factors that require special attention but have been overlooked in previous studies focused on estimating vegetation status using optical signal trajectories.Previous studies also found changes in canopy NIR reflectance caused by the exchange of older leaves for newer leaves [23,24,27].These studies considered the differences in LOPs between senesced old leaves and new leaves, but they ignored the simultaneous aging process of new leaves and mature leaves with a life span of more than three years.In this study, the observed seasonality of the canopy NIR reflectance of evergreen forests was due to age-dependent leaf optical properties at the NIR band, including both new leaf maturation and mature leaf aging, with little changes in total leaf area.New leaf maturation mainly increases NIR reflectance at the peak time in summer (May to October), and mature leaf aging mainly reduces NIR reflectance during early spring (January to June).We have already excluded differences caused by changes in the solar zenith angle [10], and our seasonally moist sub-tropical forests are free of drought impact [24,[64][65][66] or the impact of epiphylls on mature leaves in tropical forests [23,60,67], which are highly debated factors causing changes in NIR reflectance in numerous studies, and the same seasonality pattern is observed in EVI because of its strong dependence on NIR reflectance.
It is important to note that the forest canopy is composed of mixed-age leaves.Mature leaves account for a major population, while a higher proportion of young leaves are observed by the Landsat satellite because of their distribution at the top of the canopy.Both new and mature leaves influence canopy NIR reflectance, and their impacts change with new leaf expansion.New leaf expansion changes the proportions of young and mature leaves within a forest canopy over the season, thus, further changes the canopy spectral signatures.

Potential Implications for Photosynthesis
Young leaves exhibit a higher photosynthetic capacity than older leaves, and it is therefore essential to track changes in the age-structure of leaves in the canopy, which could substantially improve the modeling of the seasonal dynamics of photosynthesis.Xiaoquan and Deying [68] studied the leaf-age effects on the photosynthetic characteristics of 18-year-old Cunninghamia lanceolata stands and found similar trends in leaf photosynthetic rate at different ages: the photosynthesis rate of newly flushed leaves increases during the maturation process and then decreases at the end of the growing season, while the photosynthesis rate of mature leaves (1 and 2 a) decreases as time progresses.This study suggested that considering the temporal variation in the LOPs of new leaves matters significantly in understanding the seasonal trend of canopy spectral signatures.Ignoring variations in leaf-scale properties may mislead our interpretations of seasonality of RS signals from evergreen forests where different-age leaf cohorts coexist in the canopy.This study provides new evidence of the importance of considering the phenology of LOPs at different ages and which could contribute to more accurately modeling of photosynthesis [69].

Implications of LOPs and Canopy Structure
At the canopy scale, seasonal variations in new leaf optical properties were shown to be the dominant factor producing seasonal variations in canopy reflectance and altering NIR to red ratios, independent of changes in other canopy attributes.However, our study explicitly focused on dense evergreen forests with stable LAI (LAI > 3).The parameters that drive the GORT model may vary on different temporal scales.For example, crown size and tree height vary annually, while LAI, leaf spectra, and sun geometry vary seasonally.We can explore the signal differences caused by LAI and canopy structure by comparing the annual differences in the GORT model outputs.We only applied the annual structural parameters from 2005 to 2015 (there were no data in 2009 and from 2012 to 2014) to drive the the GORT model.In addition to annual stable structural parameters, the SZN, LAI and canopy average LOPs were updated monthly.Although the LAI seasonality and structure parameters varied from year to year, NIR reflectance showed similar seasonal trajectories with little annual difference, as illustrated by the vertical error bar in Figure 11.There may be two explanations for this limited variation.First, there is a lack of LAI seasonality or the canopy is too dense with a high LAI; thus, small changes in leaf area are not as sensitive compared with the variations in LOPs.The results might be different for young stands, especially before canopy closure.Second, canopy structural parameters do not change significantly during the new leaf maturation process; thus, canopy structures have limited impacts on canopy reflectance during this period.As illustrated in other studies [23], LAI is also an important factor contributing to the seasonality of NIR reflectance.
This study quantified the effects of leaf age on canopy reflectance in mature evergreen forests to first explore the possible impacts on leaf quality when total leaf area does not exhibit significant seasonal changes.Given the findings of this case study, future studies in deciduous forests, where seasonal variations in LAI are more significant and leaf-age groups are less complicated, would allow further validating of our findings and apply our findings to better quantify the canopy spectral signatures to understand the role of forests in terrestrial ecosystems.

Conclusions
This study evaluated the effect of LOPs for leaves at different ages on the canopy spectral signature variation during the growing season for Chinese fir stands.For closed canopy evergreen stands with relatively stable LAI, new leaves exerted disproportional influence on the canopy season spectral signature due to the spatial distribution of the new leaves in the top and outer canopy.The most direct implications of our results are related to ecological or physiological studies that utilize remote sensing, and our findings provide a promising potential to improve the interpretation of RS signals.The most significant finding of this study at the leaf scale is the increase in the NIR transmittance of mature leaves and the increase in the NIR reflectance of new leaves.To date, most of the identified contributions of leaf age to the variation in leaf optical properties have involved changes in the differences between old leaves and new leaves.Simultaneous monitoring of new and mature LOPs with season, however, has not been previously documented.In this study, we observed an approximately 11% increase in NIR reflectance (0.05 unit) for new leaves and a 35% increase in NIR transmittance for mature leaves during the growing season.
Variations of LOPs at the leaf scale have significant impacts at the canopy scale, and contribute to seasonality of canopy NIR reflectance and EVI2.Due to the complexity of forest ecosystems, analyses based on field data alone cannot provide guidance in the interpretation of RS signals.Conversely, studies based on modeling alone, without proper ground measurements of the key factors driving canopy signal variation, can be misleading.This study combined field observations with the GORT model to elucidate the effects of leaf age on canopy-scale reflectance signals.We demonstrated that, in addition to sun-sensor geometry, the effects of leaf aging on LOPs were the major factor contributing to the seasonality of canopy reflectance for the Chinese fir stands:

•
New leaf maturation is the main factor contributing to seasonality of canopy signals (NIR REF and EVI2), because of the distribution of these leaves in the top and outer canopy, as well as their increasing proportions with leaf growth.A small increase (0.05 unit) in new leaf NIR reflectance results in a significant increase in canopy NIR reflectance from spring to summer, while a decrease in new leaf NIR transmittance from August to October causes a decreasing trend in canopy NIR reflectance in autumn and winter.

•
Mature leaf aging is another factor contributing to the seasonality of the canopy signals (NIR REF and EVI2) because of the significant proportion of mature leaves in the canopy.Mature leaf NIR transmittance is greater during the growing season than off the growing season.This difference in leaf TRA causes an increased difference in canopy reflectance between winter and summer.
Thus, the effects of leaf age cannot be ignored when conducting time series analyses using RS data for the evergreen needle leaf forests.before and after year 2006.This appendix provides details on LAI data acquisition and processing to produce consistent LAI time-series.

Appendix A.3. Leaf Area Proportion and Distribution in Crown
rooftop facing toward our study sites and automatically records the above-canopy radiation every 5 min.Ninety-degree view caps were used on both units to avoid the influence of other objects on the sensors and made these measurements comparable with those from the DHP methods in four cardinal directions.The LAI-2000 measurements were taken under diffuse sky conditions in the early morning or after sunset.The TRAC method requires direct solar radiation, so we took TRAC measurements during midday at a constant walking pace along several parallel transects, which were perpendicular to the direction of tree stem shadows.Distance markers were registered every 5 m.The TRAC data were processed by TracWin software, which calculates LAI e and the clumping index (Ω).
During August 2015, we used LAI-2000 instruments to measure the LAI in the WS2 and WS3 plots.We revisited these plots in June 2016, taking measurements with both the TRAC and LAI-2000 instruments.In June 2016, we surveyed the ZH1 and FZ1 sample plots and used the LAI-2000 instruments to measure the LAI at the same locations as measured by the long-term DHP methods.In comparison to the LAI-2000 instrument, the accuracy of DHP is affected by photograph exposure settings because these settings impact the ratio of green leaves to sky.An LAI-2000 unit was operated subsequently at the same locations for comparison with DHP, and then we convert LAI e to true LAI values.
Additional correction parameters were required to convert LAI e to LAI t .The needle-to-shoot area ratio (γ e ) was measured using destructive sampling conducted in August 2015.The clumping index Ω was measured by the TRAC instrument.Additionally, the woody-to-total area ratio (α) for Cunninghamia lanceolata was derived from destructive sampling used to calculate biomass.

Appendix B.4 Unifying LAI Measurements Using Different Methods
LAI DHP was corrected using LAI e measurements from the LAI-2000 units to decrease the underestimation caused by the automatic exposure problem present in the DHP method.We made comparisons between the LAI as measured using DHP and LAI-2000 in the same locations in ZH1 and FZ1 to estimate the system bias (ε) of DHP.We applied a fixed bias instead of a regression relationship between DHP and LAI-2000 because our study site is pure Cunninghamia lanceolata plantations and the LAI values fall within a small range.Thus, it is not possible to build a robust regression relationship with limited LAI variations using the field measurements made on our study sites.Although previous studies have built and applied regression relationships between DHP and LAI-2000 LAI methods [71], we chose to estimate the system error of DHP with a fixed value in our study site.The reasons for this choice are listed as follows: first, regression relationships are usually site specific, lack universality and cannot be applied to other places.We used the existing regression relationship to our study sites, but the result was of poor quality.The regression relationship built by [71] resulted in a mean absolute error of 1.0 LAI when compared with field LAI-2000 measurements.Second, the accuracy of the regression relationship may vary in different ranges of values.For example, if the regression relationship fits well in the low-value range but fits poorly at high values, then when we apply it to other places with mainly high values, the relationship will fail.At this study site, the LAI is quite stable; as shown in Table A2, the standard deviation is quite small (approximately 0.25 as measured by LAI-2000), and the maximum and minimum LAI-2000 measurements within our study sites are quite close (3.19-3.36 for ZH1, and 3.18-3.3for FZ1).Thus, correcting for the systematic error is a sufficient and reliable method for converting LAI values measured by the DHP method to LAI measured with LAI-2000.According to Table A2, the system bias was set to 1.83 for the DHP method.The correction parameters required to convert LAI e to LAI t were set as follows: the needle-to-shoot area ratio (γ e ) was set to 1.1 according to the results of the destructive sampling conducted in August 2015.The clumping index Ω was produced by the TRAC instrument and set to 0.8.The value of the woody-to-total area ratio (α) for Cunninghamia lanceolata was derived from the destructive samples used to calculate biomass and was set to 0.2.An illustration of the data processing workflow may be useful for understanding how we unified the LAI field data measured by different instruments, as shown in Figure A5.According to Table A2, the system bias was set to 1.83 for the DHP method.The correction parameters required to convert to were set as follows: the needle-to-shoot area ratio ( ) was set to 1.1 according to the results of the destructive sampling conducted in August 2015.The clumping index was produced by the TRAC instrument and set to 0.8.The value of the woodyto-total area ratio ( ) for Cunninghamia lanceolata was derived from the destructive samples used to calculate biomass and was set to 0.2.An illustration of the data processing workflow may be useful for understanding how we unified the LAI field data measured by different instruments, as shown in Figure A5.
This study focused on two National Research Stations of Forest Ecosystems in Huitong county, Hunan province, China, as shown in Figure 1.The first station (Station 1, S1) is located at 26 • 40 N, 109 • 26 E, and the second station (Station 2, S2) is located at 26 • 50 N, 109 • 45 E. S1 was established in 1983, and S2 was established in 1996.

Figure 1 .
Figure 1.These four study plots are covered by Cunninghamia lanceolata (also known as Chinese fir) plantations, which were replanted after clear-cutting.

Figure 2 .
Figure 2. The four scene components used in the geometric optical (GO) model.C is the sunlit tree crown; T is the shaded tree crown; G is the sunlit background; and Z is the shaded background.

Figure 2 .
Figure 2. The four scene components used in the geometric optical (GO) model.C is the sunlit tree crown; T is the shaded tree crown; G is the sunlit background; and Z is the shaded background.

Figure 3 .
Figure 3. Flow chart of LOPs retrieval and application method.Figure 3. Flow chart of LOPs retrieval and application method.

Figure 3 .
Figure 3. Flow chart of LOPs retrieval and application method.Figure 3. Flow chart of LOPs retrieval and application method.

Figure 4 .
Figure 4. Results for the first (A) and second (B) global sensitivity analyses of the GORT model in the third (red) and fourth (NIR) Landsat bands for the following parameters: LAI (leaf area index), h1 (lower boundary of canopy center height), h2 (upper boundary of canopy center height),(tree stem density (trees/ha)), r (crown radius), (leaf reflectance), (leaf transmittance), and (background reflectance).

Figure 5 .
Figure 5. Single order sensitivity analysis of the impacts of sensitive parameters on GORT model outputs for the red band (A) and the NIR band (B).Only one of the model sensitive parameters was adjusted each time to study variations in the model output spectral signatures.

Figure 4 .Figure 4 .
Figure 4. Results for the first (A) and second (B) global sensitivity analyses of the GORT model in the third (red) and fourth (NIR) Landsat bands for the following parameters: LAI (leaf area index), h1 (lower boundary of canopy center height), h2 (upper boundary of canopy center height), λ (tree stem density (trees/ha)), r (crown radius), r L (leaf reflectance), t L (leaf transmittance), and r G (background reflectance).

Figure 5 .
Figure 5. Single order sensitivity analysis of the impacts of sensitive parameters on GORT model outputs for the red band (A) and the NIR band (B).Only one of the model sensitive parameters was adjusted each time to study variations in the model output spectral signatures.

Figure 5 .
Figure 5. Single order sensitivity analysis of the impacts of sensitive parameters on GORT model outputs for the red band (A) and the NIR band (B).Only one of the model sensitive parameters was adjusted each time to study variations in the model output spectral signatures.

Figure 7 .
Figure7.Field spectral curve data (mean ± 1 s.d.) in shortwave bands for new (A1-E1) and mature (A2-E2) needle samples (n) from Chinese fir collected on 5 May, 24 June, 28 July, 14 September and 13 October 2017.From May to October, n = 22, 15, 26, 30 and 39 respectively, and each sample includes 5 to 8 leaves.Reflectance data ( ) for each month are presented as the lower set of curves within each plot, while transmittance data ( ) are presented as the upper set of curves.Absorption characteristics are depicted based on the area between the upper and lower set of curves.SD values for new leaf are depicted in the upper and lower dash lines, to avoid overlap with , while SD values for all other and measurements are depicted in gray buffed areas.

Figure 7 .
Figure7.Field spectral curve data (mean ± 1 s.d.) in shortwave bands for new (A1-E1) and mature (A2-E2) needle samples (n) from Chinese fir collected on 5 May, 24 June, 28 July, 14 September and 13 October 2017.From May to October, n = 22, 15, 26, 30 and 39 respectively, and each sample includes 5 to 8 leaves.Reflectance data (r L ) for each month are presented as the lower set of curves within each plot, while transmittance data (t L ) are presented as the upper set of curves.Absorption characteristics are depicted based on the area between the upper and lower set of curves.SD values for new leaf t L are depicted in the upper and lower dash lines, to avoid overlap with r L , while SD values for all other r L and t L measurements are depicted in gray buffed areas.

Figure 8 .
Figure 8. Retrieved results for leaf reflectance (r L ) and transmittance (t L ) at the canopy scale in the NIR band (A1,A2) and RED band (B1,B2).

Figure 9 .
Figure 9. Seasonal variations in the leaf properties of the canopy, including (A) total leaf area: average LAI from stand age 22 to 33 (2005 to 2015); (B) leaf proportion: the percentage of new leaves and mature leaves in the canopy measured by Zhongkun et al. [50]; (C) weights for new leaves (w1) and mature leaves (w2) from stand age 22 to 33 (2005 to 2015).

Figure 9 .
Figure 9. Seasonal variations in the leaf properties of the canopy, including (A) total leaf area: average LAI from stand age 22 to 33 (2005 to 2015); (B) leaf proportion: the percentage of new leaves and mature leaves in the canopy measured by Zhongkun et al. [50]; (C) weights for new leaves (w1) and mature leaves (w2) from stand age 22 to 33 (2005 to 2015).

Figure 10 .
Figure 10.Validation of estimated new leaf LOPs in the NIR band (A1) and red band (A2) and the constructed mature leaf LOPs trajectories (B1,B2).

Figure 10 .
Figure 10.Validation of estimated new leaf LOPs in the NIR band (A1) and red band (A2) and the constructed mature leaf LOPs trajectories (B1,B2).

Figure 11 .
Figure 11.Aging effects of new leaves and mature leaves on seasonality of canopy signal trajectory in NIR band (A1), red band (B1) and EVI2 (C1).With the exception of the differences in leaf optical parameters, all other input parameters for the GORT model are the same in three circumstances.We evaluated the simulated canopy signature with Landsat observations (A2-C2).Both the simulated results and Landsat observations are at a monthly step with the mean value and s.d.(from 2005 to 2015).

Figure 11 .
Figure 11.Aging effects of new leaves and mature leaves on seasonality of canopy signal trajectory in NIR band (A1), red band (B1) and EVI2 (C1).With the exception of the differences in leaf optical parameters, all other input parameters for the GORT model are the same in three circumstances.We evaluated the simulated canopy signature with Landsat observations (A2-C2).Both the simulated results and Landsat observations are at a monthly step with the mean value and s.d.(from 2005 to 2015).

Figure A1 .Figure
Figure A1.Measured parameters for sample trees selected to build regression relationships for the structural parameters of the tree crown and other measures, including tree height and DBH.

Figure A1 .
Figure A1.Measured parameters for sample trees selected to build regression relationships for the structural parameters of the tree crown and other measures, including tree height and DBH.

Figure A1 .Figure
Figure A1.Measured parameters for sample trees selected to build regression relationships for the structural parameters of the tree crown and other measures, including tree height and DBH.

Figure A2 .
Figure A2.(A) Measured crown width (2R) for 40 sample trees with DBH in the north-south direction (CW NS ) and east-west direction (CW EW ); (B) Regression relationships between DBH and crown width (2R).

Figure A3 .
Figure A3.(A) Measured height values for 40 sample trees with increasing DBH.The measured data include tree height (H1), height of crown center (h) and height under crown (H2); Regression relationships between DBH and tree height (B) and between tree height and crown center height (C).

Figure A4 .
Figure A4.The vacuum side of one group of leaf leaves prepared for spectral measurement.Leaf samples were collected on 15 July 2016.

Figure A3 .
Figure A3.(A) Measured height values for 40 sample trees with increasing DBH.The measured data include tree height (H1), height of crown center (h) and height under crown (H2); Regression relationships between DBH and tree height (B) and between tree height and crown center height (C).

Figure A3 .
Figure A3.(A) Measured height values for 40 sample trees with increasing DBH.The measured data include tree height (H1), height of crown center (h) and height under crown (H2); Regression relationships between DBH and tree height (B) and between tree height and crown center height (C).

Appendix A. 3 .
Leaf Area Proportion and Distribution in Crown

Figure A4 .
Figure A4.The vacuum side of one group of leaf leaves prepared for spectral measurement.Leaf samples were collected on 15 July 2016.

Figure A5 .
Figure A5.Flowchart of LAI field data processing procedures.Before January 2007, monthly LAI was measured using a CI-110 instrument, and a DHP instrument was used afterwards.

Figure A5 .
Figure A5.Flowchart of LAI field data processing procedures.Before January 2007, monthly LAI was measured using a CI-110 instrument, and a DHP instrument was used afterwards.

Table 2 .
Prior knowledge of input parameters and corresponding value ranges used in sensitivity analysis.Adjustable parameters.* The prior knowledge is the initial value of parameters; the lower and upper limit values are the thresholds of the parameters.
4.2.Optical Properties of Individual New and Mature Leaves 4.2.1.Leaves at Different Ages

Table 2 .
Prior knowledge of input parameters and corresponding value ranges used in sensitivity analysis.Adjustable parameters.* The prior knowledge is the initial value of parameters; the lower and upper limit values are the thresholds of the parameters.
4.2.Optical Properties of Individual New and Mature Leaves4.2.1.Leaves at Different Ages Remote Sens. 2018, 10, x FOR PEER REVIEW 17 of 28 sites.Assuming that the LOPs of mature leaves do not change, if we apply mature leaf LOPs measured at the peak time in summer (NIR TRA = 0.34) to other times, we might overestimate canopy NIR REF in winter; and if we apply mature leaf LOPs measured in early spring (NIR TRA = 0.26), we might underestimate canopy NIR REF in summer.

Table A1 .
Leaf area and proportion of leaves at different ages located at different crown positions for Cunninghamia trees in Hunan.
Appendix A.4. Leaf Sample for SVC Measurements

Table A1 .
Leaf area and proportion of leaves at different ages located at different crown positions for Cunninghamia trees in Hunan.
Appendix A.4 Leaf Sample for SVC Measurements

Table A1 .
Leaf area and proportion of leaves at different ages located at different crown positions for Cunninghamia trees in Hunan.
Appendix A.4. Leaf Sample for SVC Measurements

Table A2 .
Comparison of effective LAI derived from digital hemispherical photography (DHP) and LAI-2000 methods.The LAI were measured at 5 points, in four directions, and repeated two to three times for every plot.The system bias (ε) was calculated by (LAI-2000 LAI e -DHP LAI DHP ). *