Narrowband Bio-Indicator Monitoring of Temperate Forest Carbon Fluxes in Northeastern China

Developments in hyperspectral remote sensing techniques during the last decade have enabled the use of narrowband indices to evaluate the role of forest ecosystem variables in estimating carbon (C) fluxes. In this study, narrowband bio-indicators derived from EO1 Hyperion data were investigated to determine whether they could capture the temporal variation and estimate the spatial variability of forest C fluxes derived from eddy covariance tower data. Nineteen indices were divided into four categories of optical indices: broadband, chlorophyll, red edge, and light use efficiency. Correlation tests were performed between the selected vegetation indices, gross primary production (GPP), and ecosystem respiration (Re). Among the 19 indices, five narrowband indices (Chlorophyll Index RedEdge 710, scaled photochemical reflectance index (SPRI)*enhanced vegetation index (EVI), SPRI*normalized difference vegetation index (NDVI), MCARI/OSAVI[705, 750] and the Vogelmann Index), and one broad band index (EVI) had R-squared values with a good fit for GPP and Re. The SPRI*NDVI has the highest significant coefficients of determination OPEN ACCESS Remote Sens. 2014, 6 8987 with GPP and Re (R2 = 0.86 and 0.89, p < 0.0001, respectively). SPRI*NDVI was used in atmospheric inverse modeling at regional scales for the estimation of C fluxes. We compared the GPP spatial patterns inversed from our model with corresponding results from the Vegetation Photosynthesis Model (VPM), the Boreal Ecosystems Productivity Simulator model, and MODIS MOD17A2 products. The inversed GPP spatial patterns from our model of SPRI*NDVI had good agreement with the output from the VPM model. The normalized difference nitrogen index was well correlated with measured C net ecosystem exchange. Our findings indicated that narrowband bio-indicators based on EO-1 Hyperion images could be used to predict regional C flux variations for Northeastern China’s temperate broad-leaved Korean pine forest ecosystems.


Introduction
Forest ecosystems play an important role in decreasing atmospheric CO2 concentrations and mitigating global climate change [1,2].The carbon (C) sequestration by terrestrial ecosystems and the prediction of climate change impacts on the C balance in China have been estimated using inventory and modeling methodologies during the previous two decades at the national, regional, and biome spatial scales.C sequestration in forest ecosystems has been estimated at 3.26 to 9.11 Pg C (with a mean of 5.49 Pg C) [3].In order to improve the accuracy of China's terrestrial C sequestration and storage estimates a combination of national inventory field plot measurements, point data from eddy covariance C flux tower stations, spatial scaling utilizing remote sensing techniques, and ecosystem model simulations were identified as means to fully account for the C budget.Forest ecosystem C fluxes have been measured using several methods, including closed chambers, a global network of eddy covariance (EC) techniques at the stand spatial scale [4,5], and aircraft and spaceborne remote sensing platforms and sensors.Forest C fluxes have been spatially applied using these methods from individual canopies to ecosystem and regional scales.The linkage of forest ecosystem fluxes measured with EC methods and satellite remote sensing data are some of the most promising methods for scaling up the spatially limited coverage of eddy covariance sites in China that will facilitate the development of regional estimates of C fluxes [6,7].
Remote sensing techniques provide a useful methodology to study patterns in ecosystem processes at different spatiotemporal scales [8].Methodologies proposed for regional monitoring of vegetation C fluxes rely on vegetation indices (bio-indicators), which are related to tree canopy structure and physiology.Historically, bio-indicators were derived from broadband remote sensing data including the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI).These indices were developed to avoid saturation at larger leaf area index (LAI) values and to be resistant to atmospheric effects [9].NDVI, EVI, and other multispectral vegetation indices were proposed in several studies as a proxy for LAI, canopy structure, green biomass, percent green cover, and the fraction of absorbed photosynthetically active radiation (fAPAR) [10][11][12][13].
The estimation of ecosystem C fluxes utilizing a data-driven [14] and satellite remote sensing approach was previously used to extend the net ecosystem exchange (NEE) for 1km by 1km cells across the conterminous U.S. for a 10-year period from 2000 to 2009 [15].Global terrestrial C fluxes have also been up-scaled using a machine learning technique model tree ensembles (MTE), based on climate and meteorological data, land use data, and remote sensing indices.Trained MTEs were used to develop global flux fields at a 0.5° by 0.5° spatial resolution and a monthly temporal resolution from 1982 to 2008 [16].Cross validation analyses demonstrated good performance of MTEs in predicting among site flux variability and seasonal patterns.Correlations between NEE and broadband vegetation indices (e.g., NDVI) were not as robust [17,18].These results indicate that utilizing remote sensing techniques to estimate NEE is difficult using broadband indices, and that gross primary production (GPP) is more important to inter annual variability in NEE than ecosystem respiration (Re).
The GPP from Moderate Resolution Imaging Spectrometer (MODIS) has been evaluated using regional tower EC flux observations in North America.Spring season estimates of GPP were poorly correlated to tower EC flux data owing to the relatively rapid onset of leaf-out [19].Several previous studies have shown that vegetation indices based on narrow spectral bands significantly improve the prediction of vegetation biophysical characteristics by reducing the saturation of spectral indices at high biomass values [20,21].A recent study [13] has demonstrated that chlorophyll content, light use efficiency, and chlorophyll fluorescence narrowband indices have more advantages for monitoring GPP seasonal variations than NDVI and EVI broadband indices associated with the structure of canopy in an evergreen canopy with small physiological changes.
Hyperspectral remote sensing advances have allowed the development of smaller, more stable, accurately calibrated sensors that can quantify biophysical properties on the basis of more exact spectral absorbing and scattering characteristics of the land surface [22,23].Hyperspectral vegetation indices (HVIs) take advantage of fine scale spectral sampling to develop a narrowband equivalent of a broadband VI.A potentially greater contribution of hyperspectral systems is their ability to create new indices that incorporate wavelengths not sampled by any broadband system and to quantify absorptions that are specific to important biochemical and/or biophysical quantities of vegetation [24].This technology advancement provides some opportunities for assessing additional narrowband indices, including PRI (photochemical reflectance index), NDNI (normalized difference nitrogen index), a series of related biochemical component and red edge indices (e.g., leaf chlorophyll index, chlorophyll vegetation index, canopy chlorophyll content index, and Dmax, the maximum derivative of reflectance in 650-750 nm).The hyperspectral data are primarily derived from space-borne sensors (e.g., EO-1 Hyperion) [25,26], airborne sensors (e.g., AVIRIS, Carnegie Airborne Observatory AToMS) [27,28], and ground vegetation spectral measurements collected from FluxNet EC towers [29].Although it is underutilized in the monitoring of C fluxes, hyperspectral remote sensing has much to offer toward the improved accuracy of global carbon exchange models, e.g., monitoring light use efficiency (LUE) [30].
The Hyperion Imaging Spectrometer aboard the EO-1 spacecraft [31] focuses on hyperspectral data for assessing land cover/land use, mineral resources, coastal processes, and earth/atmospheric processes that are widely applied in many aspects of ecology.Hyperion data have been used for mapping plant canopy biochemistry (e.g., nitrogen concentration) [32].It has also been employed to test indices for estimating basal area and canopy cover of savannah woodland in central Brazil [33], and to analyze the variation in chlorophyll content and LAI with bands from the red edge of the vegetation spectrum [34].
Hyperion data are also being used to monitor ecosystem C fluxes.Earth Observing System (EOS) sites provide key data for monitoring properties of major ecosystem types associated with high spectral resolution measurements (≤10 nm, 400-2500 nm) to provide synoptic evaluation of the factors significantly affecting the ability of the vegetation to sequester C and to reflect radiation due to changes in vegetation pigment and water content, structural and chemical composition.Preliminary results suggest a strong correlation between CO2 fluxes and a number of bio-indicators associated with pigment content, regardless of vegetation type and site [35].
During the previous two decades, there has been an increase in the number of global CO2 flux sites conducting continuous year-round monitoring.Scaling-up fluxes from site-specific tower measurements to regional spatial scales has been the focus of methods development for remote sensing and mathematical modeling.Measurements from EC flux towers are likely to have a C flux estimate bias when they are scaled up to regional and global assessments.Although broadband indices have been employed in large-scale studies, they do not consistently provide detailed information of spectra and spatial considerations.Meanwhile, narrowband bio-indicators need to be investigated for their utility in monitoring ecosystem features [36].In this study, we employed EO-1 Hyperion data and measured EC flux tower data in a pine and broad-leaved mixed forest to evaluate narrowband indices for monitoring C fluxes in Northeastern China.The primary goals of this study were: (i) To evaluate effectiveness of narrowband bio-indicators for monitoring temporal variations of forest C fluxes; (ii) To map regional patterns of C fluxes using the optimal bio-indicators and validate inversed GPP spatial patterns; (iii) To develop methods for the use of narrowband bio-indicators for hyperspectral remote sensing of terrestrial ecosystem C fluxes across forest stand and regional spatial scales.

Study Area
The study area was located on the northern slope of Changbai Mountain in Jilin Province, Northeast China (41°41′49″N-42°25′18″N, 127°42′55″E-128°16′48″E) (Figure 1).The annual average temperature decreases with elevation and ranges from 2.8 °C at 750 m above mean sea level (AMSL) to −0.75 °C at 1400 m AMSL.Annual precipitation increases with elevation, with an annual mean of 670 mm at 750 m to 870 mm at 1400 m AMSL Eighty percent of the total annual precipitation occurs between June and September [37].The soil in the study area is classified as dark brown forest soil originating from volcanic ashes [38].The vegetation near the tower is an approximately 200-year-old, uneven-aged, mixed forest of Korean pine (Pinus koraiensis Sieb.et Zucc.),Amur linden (Tilia amurensis Rupr.), Painted maple (Acer mono Maxim.),Manchurian ash (Fraxinus mandshurica Rupr.), Mongolian oak (Quercus mongolica Fisch.Ex Ledeb.), and 135 understory species.The study area includes additional land use types, e.g., residential, commercial, and agricultural.
The EC flux tower is located at the Chinese Academy of Sciences' Forest Ecosystem Experimental Station of the Changbai Mountains (42°24′9″N, 128°05′45″E, elevation 738 m AMSL).Mass and energy fluxes were measured continuously at the site with an EC system mounted at a height of 40 m.The system consists of a triaxial ultrasonic anemometer (CAST3, Campbell Scientific, USA) and a fast response open-path CO2/H2O infrared gas analyzer (Li-7500, LiCor Inc., Lincoln, NE, USA) mounted on an arm extending 3 m from the tower.The anemometer and the infrared gas analyzer (IRGA) provide digital outputs of the three wind components, sonic temperature, water vapor, and CO2 density.Raw data were recorded at a frequency of 10 Hz and stored every 10 s as the 30 min averages into separate files [39].

C Fluxes Data
The C flux data included three variables (GPP, Re, and NEE), but only NEE was measured directly using EC instrumentation.The GPP and Re were estimated by NEE decomposition techniques using the methodology below.
The NEE time series were decomposed into GPP and Re according to the following four equations: where NEE at no light (QPPFD = 0) is estimated as RD (the mean of daytime respiration of the time window), a1is the maximum photosynthetic uptake, and a2 is the light level (QPPFD) corresponding to half the maximum photosynthesis rate.
To estimate daytime ecosystem respiration, nighttime NEE data were calculated according to the Arrhenius type function [40] below: where R is the gas constant(8.314J•K −1 •mol −1 ), Tref is the reference temperature constant of 283.15 K, Tk is soil temperature (K), and Ea is the activation energy in J•mol −1 [38].GPP was calculated as the difference between NEE and the total ecosystem respiration (Re).
Daily ecosystem respiration Re is calculated as the sum of daytime respiration (Re,day) and nighttime respiration (Re,night): Nighttime net exchange flux is assumed to be equivalent to nighttime ecosystem respiration.Temperature response models obtained from nighttime net exchange flux (Re,night) were applied at daytime to estimate daytime ecosystem respiration (Re,day) [38].

Hyperspectral Data
Table 1 summarized the acquiring DOY (day of year) of 18 cloud free L1Gst Hyperion products.The Hyperion L1Gst products were acquired for 196 calibrated wavebands from the 242 available wavebands and excluded the zero data, VNIR/SWIR overlap, and the strong water vapor absorption bands 1-7, 58-78, and 225-242.Bands 8-57 and 79-224 included in the analysis are in the 427-925 nm (visible and near-infrared (VNIR)) and 933-2396 nm (shortwave infrared (SWIR)) regions of the electromagnetic spectrum, respectively.The signal-to-noise ratio (SNR) for Hyperion changes from 190 to 40 as the wavelength increases.Wavelength windows from 1356 to 1466 nm and 1810 to 2012 nm were eliminated due to strong disturbance by atmospheric vapor and high noise [41,42].The main spectral regions chosen in our study were 490-1336 nm, 1477-1790 nm, and 1981-2355 nm.Hyperion Level 1 data products were calibrated to sensor radiance (Wm Radiometric corrections were done for Level 1 radiance images to achieve a uniform calibration gain for each band and enhance interpretation of radiance spectra.Geometric corrections were also conducted for Level 1 L1Gst products.Bands were transformed to canopy reflectance images using the MODTRAN4-based atmospheric correction algorithm Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) built into the Environment for Visualizing Images (ENVI) software.
Reflectance spectra were then extracted and averaged from the 1 km 2 square region, which included the eddy covariance tower at its center.The spectral data were extracted and averaged for the same region.

Data Analysis
For our study, data were acquired from the CBS EC tower and the EO-1 Hyperion sensor.A semi-empirical statistical modeling approach was used to correlate hyperspectral remotely sensed surface characteristics indicated by bio-indicators to measure forest C fluxes data across different satellite overpass times.Partial least square regression (PLSR) analysis was used to explore sensitive bands and indices of observed GPP and Re.Ordinary linear regression (OLR) fitting was employed to build models.The optimal model was selected by comparing their correlation coefficients to GPP and Re.Using the optimal model, we then derived the regional maps of forest C fluxes that depict the daily C flux values for the temperate broad-leaved Korean pine forest landscapes.We validated inversed GPP by comparing our regional GPP map to the simulated regional GPP with the VPM and BEPS models, and MODIS MOD17A2 GPP products in different seasons, and explored the potential of utilizing narrowband bio-indicators.The detailed methods of key processing are described below.The study flowchart is presented in Figure 2.

Spatial Agreement between the Fluxes Tower Footprint and the Image Region
The footprint may be interpreted in analogy to the "field of view" of the instrument.The footprint was influenced by the atmospheric stability (e.g., wind speed and direction), observation height, and the underlying surface roughness.It was calculated by analytical models of C flux contributions, based on approximation to the actual physical process [43,44].The CBS flux tower footprint was defined using methods developed from a previous study [44].This study [44] reported that the footprint (distance from flux tower) of CBS Flux station ranged between 181 m to 3070 m.The footprint was delimited as a 1 km 2 square region, which included the eddy covariance tower at its center.The reflectance spectra were then extracted and averaged from the region to approximate spatial coverage agreement between the flux tower footprint and the image region.

Identification of Sensitive Wavebands and Indices to C Fluxes
A PLSR statistical method based on Principal Component Analysis (PCA) was used to determine sensitive wavebands of Hyperion images to estimate C fluxes.PLSR can analyze data with strong colinearity, noise, and numerous x-variables, and can also simultaneously process several response variables, Y [45].In the PLSR models, the importance of every x-variable both to explain x and to correlate to y was calculated in the PCA as the VIP (Variable Importance for a PLS Projection) value.VIP values larger than 1 indicate "important" x-variables, and values lower than 0.5 indicate "unimportant" x-variables.The interval between 1 and 0.5 is a gray-zone, where the importance level depends on the size of the data set.Thus, the waveband's important to GPP and Re can be considered as a selection criterion for remote sensing bio-indicators.
Similarly, sensitive bio-indicators also were selected by analyzing the VIP value for GPP and Re.The values were calculated for each xk by summing the squares of the PLS-weights, wak, weighted by the amount of y-explained in each model component, a.

The Principles of Indices Selection
Plants assimilate atmospheric CO2 through photosynthesis.The photosynthetic capacity is limited by environment factors and canopy foliar biochemical component, especially foliar nitrogen and chlorophyll content [46,47].Thus, indices related to foliar nitrogen and chlorophyll content were given priority to monitor the variation of C fluxes.
Secondly, red-edge indices have been employed to monitor vegetation characteristics [48,49].They may reflect the physical condition of the forest canopy, which can be reflected by red edge indices, can influence its photosynthetic capacity and result in C fluxes variation.Red edge indices were calculated form Hyperion data in this study.
Finally, LUE model methodology for quantifying vegetation productivity (or the net photosynthetic C gain of vegetation) can be limited by the LUE and the fraction of photosynthetically active radiation absorbed by vegetation (fAPAR) [50][51][52].We assumed that the C flux variables were a function of LUE and fAPAR.In our study, C fluxes were estimated as a function of indices, which included LUE and fAPAR.

Reflectance Curves Characteristics in Different Seasons
Vegetation reflectance varies in different wavebands and seasons.The reflectance spectrum in the growing season (Figures 3a,b,d,e) illustrates a peak in the green region and absorption features in blue and red regions.The red edge is positioned in the range of 700-720 nm wavelengths, which represents a large jump from the bottom of chlorophyll absorption (660-670 nm) to the near infrared plateau in the growing seasonal spectrum.As a result of the vegetation's high reflectance in the near infrared region and three moisture absorption valleys at 1450 nm, 1950 nm, and 2550 nm, there are two reflection peaks in the regions of 1600-1650 nm and 2150-2200 nm.
In contrast, the reflectance spectrum in the dormant season (Figure 3c,f) shows dissimilar reflectance amplitude and spectral features to the spectrum in growing season.The absorption features in the blue and red regions and reflection peaks in the green region almost disappear.The red edge features recede, the NIR shoulder decreases, and the reflection peaks flatten in SWIR region.
The NIR reflectance varies with different seasons.Some reflectance values were much higher than 0.4 during the growing season, but some were lower than 0.2 in the dormant season.

Reflectance First-Order Derivative Curves Characteristics in Different Seasons
First-order derivatives of reflectance (FDR) illustrate the variability of wavelengths in the slope of the original reflectance spectrum.The FDR curves are illustrated in the spectrum of 350-950 nm in Figure 4.The first-order derivative spectrums in different seasons were not similar, especially in the red edge region.During the growing season, FDR peaks existed in the red edge region.In different DOYs, the peak value discrepancies illustrate an important feature of peaks in the visible and the red edge regions.The first derivative of the red edge region was much higher in the growing season vs. the dormant season.

Reasons for the Spectrum Seasonal Variation
The vegetation canopy reflectance spectrum was subject to many controlling factors, notably, concentrations of canopy pigments, effects of variations in number of leaf layers (LAI), orientation of leaves (leaf angle distribution (LAD)), canopy ground coverage, presence of non-leaf elements, and areas of shadow and soil/litter surface reflectance.These variables influenced the reflectance spectrum in different wavelength regions [48].The pigment played a dominant role in visible and red edge regions (400-800 nm).Canopy structural factors have dominant effects in the near-infrared region (750-1500 nm) [53].
Due to the influence of the snow cover below and above the forest canopy, the spectrums of remote sensing pixels were mixed by the reflectance of forest canopy and snow.Snow has a high reflectance in the visible (VIS) and near-infrared (NIR) wavelengths when compared to other natural targets [54].Variations in snow reflectance can induce large changes in the scene reflectance, even though visible snow-covered ground accounted for only half of the observed area as a result of the blocking effect of the trees [55].The ratio of the snow cover to canopy cover can influence on spectrums of pixels in C flux footprint.A few of reflectance curves were much higher than the whole spectrums curves during the winter.
The existence of peaks, including the red edge peak, in the FDR curves can be attributed to the characteristic absorption spectrum of chlorophyll [56], with a possible contribution from chlorophyll fluorescence [57].The red edge peak feature can be used to identify bio-indicators.

Indices Selection
Results of the VIP values from PLSR illustrate that the most sensitive wavebands correspond to GPP and Re (Figure 5).They are located in the NIR region (740-1350 nm) and the 620-690 nm and 1980-1990 nm wavelengths, with a VIP value > 1.The most insensitive wavebands were located at 710-730 nm, 1470-1520 nm, 1740-1790 nm, and 2100-2180 nm wavelengths, with VIP value < 0.5.Other wavebands showed moderate sensitivity to GPP and Re.Based on the principles of indices and sensitive wavelength ranges selection methodologies, 11 narrowband HVIs were used for indicating C fluxes.Most of the 11 narrowband HVIs were formed by combined sensitive and insensitive bands.Eight broadband VIs were also selected for identifying the effects between broadband VIs and HVIs.VIs categories have been previously reported in the literature [58].However, we classified the 19 indices listed in Table 2 [9,35,36,[59][60][61][62][63][64][65][66][67][68][69] based on our study goals.All VIs were calculated based on corresponding Hyperion wavebands.The SPRI can be calculated from PRI by equation "scaled photochemical reflectance index (SPRI) = (PRI + 1)/2".

Seasonal Variation of All Indices
The variations of C fluxes from 2008 to 2009 are illustrated in Figure 6.The C fluxes are presented as a significant single-peaked variation characteristic in one year.C fluxes can increase to the peak values late in the growing season (summer) and decrease at the end of the growing season each year (Figure 6).The data presented in Figure 7  In order to investigate the seasonal variation of different indices, all indices were divided into different index types (Table 2) and subsets (Figure 7a-f).For different index types, seasonal change trends were variable.The indices in the same subset had a similar seasonal change trend, e.g., C fluxes (single-peaked characteristic annually), especially in the narrowband subsets of Figure 7c-f   Abbreviations: CCCI (Canopy Chlorophyll Content Index); CVI (Chlorophyll vegetation index); LCI (Leaf Chlorophyll Index); CI (Chlorophyll Index RedEdge 710).Correlation coefficients were calculated between C fluxes and the 19 indices (Table 4).These indices had various sensitivities to GPP and Re.Thirteen indices had significant correlations (p < 0.01 or p < 0.05) with GPP and Re simultaneously.The coefficients were generally higher for Re than for GPP.Most narrowband indices showed higher sensitivity to GPP and Re than most broadband indices.MCARI/OSAVI[705,750] had the highest correlation coefficients to GPP (R = 0.904) and the Vogelmann Index had the highest correlation coefficients to Re (R = 0.937).

Model Fitting and Optimal Model Selection
Based on the comprehensive comparison of large to small VIP values, six bio-indicators were selected to build and calibrate the empirical models, i.e., MCARI/OSAVI[705,750], Vogelmann Index, Chlorophyll Index RedEdge 710, EVI, SPRI*EVI, and SPRI*NDVI (Table 3).
In order to improve the coefficients of determination between the six bio-indicators and GPP and Re, linear models, polynomial quadratic, exponential growth, and logarithmic model fittings were employed to conduct empirical calibration and to develop inversion models (Figure 8).SPRI*NDVI had the highest coefficients of determination (R 2 ) to GPP and Re (R 2 = 0.86 and 0.89, p < 0.0001, respectively).
The two calibration functions listed below were used to inverse the regional forest ecosystem GPP and Re.

Inversed Spatially Distributed C Fluxes
One cloud-free Hyperion image (acquired in 17 June 2008, DOY: 166) was selected to inverse the spatial patterns of GPP, Re and NEE.The calibration functions (Equations ( 5) and ( 6)) were used to inverse GPP and Re of regions covered by vegetation.The regional NEE was then calculated by the difference of inversed GPP and Re.
GPP, Re, and NEE had similar spatial patterns, all decreasing from northeast to southwest (Figure 9).Forestland cover in the study site is a C sink in the northeast, but a C source in the southwest.This phenomenon may be correlated with the region's elevation gradient.The altitude across the study site influences local air temperature and ecosystem productivity.Across the study area, the altitude varies from 700 m AMSL to above 2000 m AMSL [37].The C fluxes were low in the land use areas with reduced vegetation (e.g., villages, rivers, and roads).

Validation of Inversion GPP
Due to the lack of spatially explicit measured data, the regional inversion GPP was validated by comparing the output of three independent models, the Vegetation Photosynthesis Model (VPM) [70], the Boreal Ecosystems Productivity Simulator (BEPS) model [71], and the MODIS MOD17A2 GPP products [72].The VPM model used two improved vegetation indices (EVI and Land Surface Water Index (LSWI)) and is capable of tracking seasonal dynamics and inter-annual variations in GPP of evergreen needle leaf forest at a sub-monthly temporal resolution.The BEPS model is a daily process model that can estimate NPP (also GPP) for large areas by using remote sensing, gridded daily meteorological data, and soil data.The algorithm of MODIS GPP (MOD17A2) is an application of the radiation conversion efficiency concept that predicts daily GPP, using satellite-derived FPAR (from MOD15) and independent estimates of PAR and other surface meteorological fields.
The GPP data output from the models have different spatial resolutions.In our study, data were reported with a spatial resolution of 30 m.The model outputs from the VPM have a spatial resolution of 500 m, and the BEPS and MOD17A2 have a spatial resolution of 1000 m.Due to the spatial resolution differences, and the uncertainty of the spatial re-sampling method, data spatial matching was problematic.We chose to contrast their spatial patterns with a combining pixel value frequency distribution analysis.One cloud-free image acquired on DOY 166 in 2008 was selected to validate our result with the three independent models.

Comparison of GPP Spatial Patterns
Figure 10 illustrates the agreement between the GPP spatial patterns and the detail features inversed by our model to the pattern inversed by the VPM model.The daily GPP significantly decreased from northeast to southwest across the rise in elevation.The spatial resolution of VPM was 500 m, and there were mixed pixels adjacent to the no vegetation regions, which may have influenced inversed GPP values in these pixels.The detail features were most consistent with our inversed results adjacent to the mask regions.Due to the coarse 1000 m spatial resolution of two other models, there were larger numbers of mixed pixels that increased uncertainty.The spatial details disappeared in the GPP patterns of BEPS and MODIS.

Validation of GPP Frequency Distribution
Although the spatial resolution and the magnitude of pixel number were variable in the four maps illustrated in Figure 10a-d, these differences did not influence the statistical histogram comparisons.The GPP frequency distribution of our inversion had double peaks, with a high peak at 7.5 gC•m −2 •day −1 and a relative low peak at 2.5 gC•m −2 •day −1 (Figure 11).There were also double peak features in the VPM and MOD17A2 outputs.No obvious double peaks were found in the BEPS output.The double peaks of the VPM result were located at 8 and 2 gC•m −2 •day −1 , which agreed with our inversion.We conclude that our model can describe the GPP spatial distribution more accurately than the other three models evaluated in our study.
Within the range of regional inversed GPP, there were differences in the outputs of the four models.The GPP of VPM was distributed in the range of 0 to 15 gC•m −2 •day −1 , the BEPS model value was distributed in the range of 2 to 17 gC•m −2 •day −1 , and the MOD17A2 had GPP values ranged from 2 to 8 gC•m −2 •day −1 .The GPP variation range in our inversed results ranged from 2 to 10 gC•m −2 •day −1 .The analysis of our model's spatial patterns combined with the pixel value frequency distribution showed promise for inversing regional GPP features.

Correlation of Bio-Indicators with NEE
The bio-indicators have significant correlations with GPP and Re, however we did not find significant correlation relationships between all bio-indicators and NEE.The non-significant correlations were a result of two factors.NEE showed a large variability during short temporal periods of several days in contrast to GPP and Re (Figure 6).Bio-indicators can only evaluate the condition of the forest canopy during the time periods of the satellite overpasses.This means that bio-indicators cannot be used to accurately evaluate the daily variation of NEE.A second reason for the observed significant correlations was that in the winter season the snow below the forest canopy has an influence on the canopy reflectance in VNIR.This resulted in low correlations of bio-indicators for estimating NEE in the winter.
To overcome these sources of error, we removed the Hyperion images acquired in winter from our model validation.Ten growing season images were used to analyze the correlation between bioindicators and NEE (Table 1).The average NEE of 5 days (NEE5) was used to determine the correlations.Eight bio-indicators had significant correlations with NEE5 (Table 5).Among the seven bio-indicators, NDNI correlated with NEE5 (R = −0.643,P = 0.045) and had a significant correlation with NEE (R = −0.805,P = 0.0050) (Figure 12).This result suggests that narrowband vegetation index with additional SWIR information might be useful to monitor the variation of forest ecosystem NEE since canopy nitrogen has been shown to promote NEE.

Performance of Narrowband Bio-Indicators
Our findings have demonstrated that the narrowband bio-indicators related to chlorophyll, red-edge, and the LUE model can be used to estimate forest productivity and respiration.In general, they have higher correlation coefficients than broadband indices.SPRI*NDVI was the best bio-indicator employed to inverse GPP and Re at the regional scale.The MCARI/OSAVI[705, 750], Vogelmann Index, Chlorophyll Index RedEdge 710, and SPRI*EVI also had high correlation coefficients with GPP and Re.
The LUE model indices exhibited the best correlation coefficients to facilitate spatially extrapolation of regional scale estimates of C fluxes.Due to the differences of NDVI and EVI in decreasing the LAI saturation effect and atmosphere influence, SPRI*NDVI may be more suitable for monitoring the vegetation productivity with relative low LAI and low atmosphere influence environments.SPRI*EVI may have more potential for high LAI vegetation types and atmosphere influence environments.SPRI*NDVI and SPRI*EVI should be tested in more detail in different ecosystems and environments in future studies.
MCARI/OSAVI[705, 750] have been proved to have good linearity with chlorophyll content and low sensitivity to non-photosynthetic element effects, mainly at low chlorophyll concentrations and low LAI.The Vogelmann Index was formed by a GPP sensitive band at 740 nm and an insensitive band at 720 nm.The ratio of the two bands can indicate the GPP and Re variations.The Chlorophyll Index RedEdge 710 was also calculated by two relevant bands in the red edge region.All three narrowband bioindicators showed good performance.

Uncertainties, Errors, and Accuracies
To validate the selection of the six bio-indicators, MCARI/OSAVI[705,750], Vogelmann Index, Chlorophyll Index RedEdge 710, EVI, SPRI*EVI, and SPRI*NDVI, and to assess the accuracy of the regional inversion GPP from our model, we compared the output of three independent models, the VPM [70], the BEPS model [71], and MOD17A2 GPP products [72].The GPP spatial pattern agreement between our model and the three independent models is illustrated in Figure 10.Statistical histogram comparisons showed good agreement between VPM and MOD17A2.The regional inversed GPP outputs of the four models were distributed within the range of 0 to 15 gC•m −2 •day −1 for VPM, 2 to 17 gC•m −2 •day −1 for BEPS, 2 to 8 gC•m −2 •day −1 for MOD17A2, and 2 to 10 gC•m −2 •day −1 for our model.Although all four models are within a similar regional inverted GPP range, improvements in model agreement were constrained by spatial resolution and the magnitude of pixel number among the four models.
Additional uncertainty and errors of model output can be attributed to the influence of snow cover on the reflectance spectrum of remote sensing pixels and the delimitation of C fluxes footprint.The ratio of the snow cover to canopy cover influences the spectrums of pixels in the VIS and NIR wavelengths.The footprint can be influenced by the atmospheric stability of wind speed and direction.Due to the variation of daily atmospheric conditions, the footprints vary based on satellite imaging times, and the fixed spatial region of extracting image spectrums vary for each model footprint.Also, the impacts of land-cover/land-use (LCLU) change on vegetation growth affected by various human activities, including residential and commercial development, reforestation, and agriculture, may introduce error in the correlation of hyperspectral remotely sensed surface characteristics indicated by bio-indicators to measure regional forest C fluxes data.

Limitations and Suggestions
There are limitations when the narrowband bio-indicators and empirical calibration models are used for routine mapping of ecosystem C fluxes.The primary limitation derives from the semi-empirical feature of this method, which requires the availability of EC tower data for calibration of narrowband bio-indicators.Semi-empirical models may lack robustness when applied to different biomes and need to be validated.Hyperspectral data are still not available in all regions where the EC towers are operational.
The other limitation is the lack of accurate validation approaches.Although the spatial pattern inversed in the study was consistent with the results of other models, the maps of C fluxes still could not be validated quantitatively due to the lack of observational data at spatial and temporal scales.
Despite these limitations, the strong correlations between bio-indicators and EC flux tower data indicate that narrowband bio-indicators are a valid approach for mapping C fluxes over large spatial scales.This approach will be improved following the completion of a new global spectral network, SPECNET, being built in China.The SPECNET will be composed of multi-angles spectrometer and eddy covariance systems.It will link hyperspectral spectrum measurements with mass and energy flux sampling across biome spatial scales.This network will enhance the linkage of EC tower fluxes and hyperspectral data across different terrestrial ecosystems and spatial scales, improving the application of narrowband bio-indicators in China.

Conclusions
In this study, a series of narrowband bio-indicators derived from EO-1 Hyperion were employed to monitor C fluxes in a Korean pine and broad-leaved mixed forest landscape on the northern slope of Changbai Mountain in northeast China for the first time.Three categories of bio-indicators, red edge indices, chlorophyll content indices, and indices based on LUE principles, demonstrated improvements for capturing the variation of GPP and Re.The scaled photochemical reflectance index (SPRI)*NDVI had the best goodness of fit for GPP and Re.The bio-indicators used for inversing regional C fluxes were SPRI*NDVI, a function of NDVI (related to canopy fAPAR) and narrowband SPRI (a dynamic efficiency factor related to canopy light use efficiency).We found a significant relationship between SPRI*NDVI and GPP and Re, R 2 = 0.86 and 0.89, p < 0.0001, respectively.The MCARI/OSAVI[705, 750], Vogelmann Index, Chlorophyll Index RedEdge 710, and SPRI*EVI had potential for utilization as bio-indicators for monitoring the variation of regional ecosystem productivity.The inversed regional average GPP was 6.22 gC•m −2 •day −1 ; the average Re was 5.78 gC•m −2 •day −1 ; and the average NEE was 0.44 gC•m −2 •day −1 on DOY 166 in 2008.These values indicate that the study area is a carbon sink region in the growing season.The inversed daily GPP spatial patterns were consistent with corresponding results from the VPM Model.
The NDNI had a significant correlation with the NEE, with R = 0.805, p < 0.01.This correlation has not been previously described in the literature.The relationships between other narrowband bio-indicators and average NEE5 were significant with p < 0.05.We believe the methodology needs additional validation with additional forest land cover types, however, the direct correlation between NEE and narrowband bio-indicators described in our study provides a simple and effective approach to monitor C fluxes variation and map the spatial patterns of regional NEE.Some limitations still exist when bio-indicators and models are used for routine mapping of C fluxes in different vegetable types.We predict that the method will become more useful with the accumulation of matched spectral and C flux data.For example, with the wide application of MODIS, SPRI*NDVI could be evaluated at large spatial scales.We believe that, in the near future, more narrowband bioindicators could be identified as an alternative method for monitoring and modeling of vegetation ecosystems at large spatial scales.

Figure 1 .
Figure 1.Location of the flux tower study area and the research station on the northern slope of Changbai Mountain in Jilin Province in Northeast China.(a) shows a China province map; (b) shows a Hyperion true color image with a non-forest land cover mask (white areas); (c) is an image of forest landscape adjacent to the CBS EC flux tower during the growing season.

Figure 2 .
Figure 2. The study flowchart.(White rectangles illustrate process modules, grey rectangles illustrate multi-source data integration modules, and ellipses illustrate the start and the result modules).

Figure 3 .
Figure 3. Hyperion average reflectance spectrum of the Korean pine and broad-leaved mixed forest in 1 km 2 square region in different DOY of 2008 and 2009.Wavelength regions where significant plant absorptions occur are indicated.(a-f) show the reflectance spectrum in different seasons in 2008 and 2009.The spectrum sample size of Figure 3a-f is 1089 pixels.

Figure 4 .
Figure 4. Hyperion average reflectance first-order derivative spectrum of the Korean pine and broad-leaved mixed forest in 1 km 2 square region in different DOY of 2008 and 2009.Wavelength regions where significant plant absorptions occur are indicated.(a-f) show the FDR spectrum in different seasons in 2008 and 2009.The spectrum sample size of Figure 4a-f is 1089 pixels.

Figure 5 .
Figure 5.The Variable Importance for a PLS Projection (VIP) values curves to gross primary production (GPP) and ecosystem respiration (Re) in the whole spectrum region.
illustrates the temporal changes of the 19 indices from 2008 to 2009 and contrasts their C flux change trends.
. NDVI, EVI, and some narrowband chlorophyll and nitrogen content indices had a single-peaked seasonal variation each year.All the red-edge indices and the LUE model indices had consistent temporal change trends with GPP and Re.The single band reflectance index did not capture the seasonal variation of C fluxes.

Figure 6 .
Figure 6.Temporal changes of C fluxes from 2008 to 2009.The red curve in the left plot illustrates the moving average line within 10 values.(a-c) show the annual variation of NEE, GPP and Re respectively from 2008 to 2009.

Figure 7 .
Figure 7. Seasonal variation (2008-2009) of 19 indices in five index categories.(a) includes four single-band indices; (b) includes four broadband ratio indices; (c,d) include five chlorophyll and nitrogen content indices; (e) includes four red-edge indices; (f) includes two light use efficiency (LUE) model indices.

4. 3 . 1 .
Importance of Bio-Indicators for C Fluxes All indices that are of importance to C fluxes (GPP and Re) are identified in VIP values in Table 3. Eleven indices have lager VIP values than 1 and are considered as effective bio-indicator indices.

Figure 8 .
Figure 8. Scatter plots of 6 bio-indicators and measured GPP and Re.All fitting models are significant at the 0.01 level (2-tailed).

Figure 9 .
Figure 9. C fluxes spatial patterns in DOY of 166 in 2008 inversed by Hyperion image with narrowband bio-indicator scaled photochemical reflectance index (SPRI)* normalized difference vegetation index (NDVI).

Figure 10 .
Figure 10.Comparison of GPP spatial patterns inversed from different models in the same time period.(a-d) illustrate the inversed results of our model, vegetation photosynthesis model (VPM), boreal ecosystems productivity simulator (BEPS), and MOD17A2, respectively.The no vegetation or cloud-covered regions are masked by white pixels in Figure 10a.

Table 1 .
The acquisition day of year (DOY) of Hyperion images used for this study.

Table 2 .
Nineteen optical indices used in the study.

Table 3 .
Importance (VIP values) of 19 indices to GPP and Re calculated from partial least square regression (PLSR).

Table 4 .
Correlations of 19 indices to GPP, Re and net ecosystem exchange (NEE).

Table 5 .
Correlation coefficients between the seven bio-indicators and NEE/NEE5 for a 200-year-old Korean pine and broad-leaved mixed forest.