Modeling Gross Primary Production of a Typical Coastal Wetland in China Using MODIS Time Series and CO 2 Eddy Flux Tower Data

How to effectively combine remote sensing data with the eddy covariance (EC) technique to accurately quantify gross primary production (GPP) in coastal wetlands has been a challenge and is also important and necessary for carbon (C) budgets assessment and climate change studies at larger scales. In this study, a satellite-based Vegetation Photosynthesis Model (VPM) combined with EC measurement and Moderate Resolution Imaging Spectroradiometer (MODIS) data was used to evaluate the phenological characteristics and the biophysical performance of MODIS-based vegetation indices (VIs) and the feasibility of the model for simulating GPP of coastal wetland ecosystems. The results showed that greenness-related and water-related VIs can better identify the green-up and the senescence phases of coastal wetland vegetation, corresponds well with the C uptake period and the phenological patterns that were delineated by GPP from EC tower (GPPEC). Temperature can explain most of the seasonal variation in VIs and GPPEC fluxes. Both enhanced vegetation index (EVI) and water-sensitive land surface water index (LSWI) have a higher predictive power for simulating GPP in this coastal wetland. The comparisons between modeled GPP (GPPVPM) and GPPEC indicated that VPM model can commendably simulate the trajectories of the seasonal dynamics of GPPEC fluxes in terms of patterns and magnitudes, explaining about 85% of GPPEC changes over the study years (p < 0.0001). The results also demonstrate the potential of satellite-driven VPM model for modeling C uptake at large spatial and temporal scales in coastal wetlands, which can provide valuable production data for the assessment of global wetland C sink/source.


Introduction
As a staggered transition zone of terrestrial ecosystem and marine ecosystem, coastal wetland constitutes chains of marshes and swamps, and provides valuable ecosystem services and benefits to human [1][2][3][4].It plays an important role in the principal interchanges of material and energy and the global sequestration of carbon (C) [1,5,6].As the key component of the C cycle, gross primary production (GPP) represents the capacity of vegetation to uptake C from the atmosphere [7][8][9][10].However, global-scale and measurement-based estimation of the historical growth in GPP is still lacking.The dynamics of GPP in coastal wetlands and its response to global climate change is poorly quantified and understood [5,11,12].Therefore, accurate quantification of GPP for coastal wetlands at various spatial and temporal scales is crucial and necessary for global C budgets assessment and climate change studies [13].
The eddy covariance (EC) method and remote sensing (RS) technology have greatly increased the opportunities for the quantification of GPP at regional or global scales.The EC method can continuously and reliably measure net exchanges of C, water and energy between atmosphere, and terrestrial ecosystems at the diurnal, seasonal, and annual scales [14][15][16][17].It can be used for calibrating and validating ecological models by providing important information that is associated with photosynthetic period and GPP at the ecosystem scale [18][19][20].However, C flux data provided by the EC method is limited by footprints (typically ranging from hundreds of meters to 1 km) and short time periods [21][22][23][24].Satellite RS technology can alternatively provide rapid, continuous, systematic, and repetitive measurement of vegetation structure and function characterization over large areas [23,25].It is an effective and powerful tool for monitoring properties of land surface and obtaining spatiotemporal information, and then estimating regional GPP and net primary production (NPP), which can be consistent with the EC tower footprint [26,27].However, there is still lacking interrelated links and bridges between EC and RS.
The satellite-based model is the interrelated links and bridges between EC measurement and RS data at regional scale.Among many RS models, the production efficiency models (PEM) are widely used in the GPP estimation study because of its lesser number of driving parameters, accurate estimation, and the easier coupling of EC and RS [23,[28][29][30][31].These models generally predict GPP by light use efficiency (ε g ), photosynthetically active radiation (PAR), the fraction of the absorbed PAR by the vegetation canopy (FAPAR), and a greenness-related Normalized Difference Vegetation Index (NDVI) [32].Recently, Xiao et al., has developed a RS-based vegetation photosynthesis model (VPM) to simulate GPP [33][34][35].VPM is a well-established GPP model by the conceptual partitioning of chlorophyll and non-photosynthetically active vegetation within the canopy to simulate GPP [36].The model used additional spectral bands, such as blue and shortwave infrared bands to calculate the Enhanced Vegetation Index (EVI) [37] and Land Surface Water Index (LSWI) [33,36,38].For example, Li et al., used the improved EVI and LSWI as input data for the VPM model to better characterize the dynamics of vegetation.They also simulated GPP more accurately than other PEM models, which only used NDVI in Qinghai-Tibetan Plateau, China [23].Recently, this model has been further developed and applied to several ecosystems, including forest, grassland, agriculture, and alpine wetland ecosystems [10,13,15,23,33,36,39,40], which demonstrated the potential to scale up in situ observations of GPP from the CO 2 flux tower sites.
As the largest wetland ecosystem in the warm temperate zone of China, the Yellow River Delta Wetland is the most active area of atmosphere-land-ocean interaction in the world [41][42][43][44].Therefore, the accurate estimation of the C source/sink function and its influence mechanism for the targeted coastal wetland will help to quantify the interactions and feedbacks between global climate change and coastal wetland ecosystems.In this study, we have combined the analysis of EC flux data with remote sensing images in the Yellow River Delta during 2009-2010.We aim to: (i) characterize seasonal dynamics of reed by analyzing the Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation indices (VIs) and CO 2 fluxes data; (ii) examine the biophysical expression of VIs related to seasonal dynamics of observed GPP data; and, (iii) evaluate the dependability and the applicability of the VPM model for simulating GPP of coastal wetland ecosystems.We also conducted a variety of uncertainty analysis for the model simulations.This study could help to deepen our understanding of Remote Sens. 2018, 10, 708 3 of 20 C budgets of coastal wetland and further improve the potential of RS technology for monitoring and simulating the dynamics of regional vegetation and C fluxes.

Description of the Study Region
Our study was carried out in the Yellow River Delta Nature Reserve (37 • 45 59" N, 119 • 09 05" E, −4 m above sea level; Figure 1).The nature reserve covers an area of 15.3 × 10 4 ha.Our EC tower was established at the coastal zone of the Yellow River Delta Nature Reserve.In this region, the dominant wind direction is the north and northwest wind.The average annual wind speed is about 4 m/s.The annual mean air temperature is 12.1 • C [41,45].The annual mean precipitation is 520 mm, which mainly concentrates in summer.The annual evaporation is 1962 mm and frost-free period is 196 days.The dominated soil type is coastal saline soil.The pH of the soil is 7.6-8.5.The dominant vegetation includes both the herbaceous plants (Phragmites australis and Suaeda salsa) and the shrub (Tamarix chinensis).The height of plants ranges from 120 to 200 cm.C budgets of coastal wetland and further improve the potential of RS technology for monitoring and simulating the dynamics of regional vegetation and C fluxes.

Description of the Study Region
Our study was carried out in the Yellow River Delta Nature Reserve (37°45′59″N, 119°09′05″E, −4 m above sea level; Figure 1).The nature reserve covers an area of 15.3 × 10 4 ha.Our EC tower was established at the coastal zone of the Yellow River Delta Nature Reserve.In this region, the dominant wind direction is the north and northwest wind.The average annual wind speed is about 4 m/s.The annual mean air temperature is 12.1 °C [41,45].The annual mean precipitation is 520 mm, which mainly concentrates in summer.The annual evaporation is 1962 mm and frost-free period is 196 days.The dominated soil type is coastal saline soil.The pH of the soil is 7.6-8.5.The dominant vegetation includes both the herbaceous plants (Phragmites australis and Suaeda salsa) and the shrub (Tamarix chinensis).The height of plants ranges from 120 to 200 cm.

CO2 Flux and Climate Data from the Eddy Flux Tower Site
The EC instrument was established at the study site since 2008 (Figure 1), with the observation of CO2 and water fluxes, wind speed, and wind direction.The EC tower equipped with an open path infrared CO2/H2O gas analyzer (Li 7500, LI-COR Inc., Lincoln, NE, USA) and ultrasonic anemometer (CSAT3, Campbell Scientific Inc., Jackson, MS, USA), which were installed 4.5 m above ground level.The sonic anemometer and the infrared CO2/H2O gas analyzer output were recorded at half-hour intervals by datalogger (CR5000, Campbell Scientific, Inc., Jackson, MS, USA), with high frequency of 10 Hz.Carbon and energy fluxes were calculated by using the half-hour covariance of vertical wind velocity and virtual temperature, water vapor density, and CO2 density [15].
We also measured some environmental factors that were used for gap filling calculations nearby the eddy tower.Photosynthetically active radiation (PAR) was observed by LI-190SB sensors (LI-COR Inc., Lincoln, NE, USA).Relative humidity and air temperature were observed by

CO 2 Flux and Climate Data from the Eddy Flux Tower Site
The EC instrument was established at the study site since 2008 (Figure 1), with the observation of CO 2 and water fluxes, wind speed, and wind direction.The EC tower equipped with an open path infrared CO 2 /H 2 O gas analyzer (Li 7500, LI-COR Inc., Lincoln, NE, USA) and ultrasonic anemometer (CSAT3, Campbell Scientific Inc., Jackson, MS, USA), which were installed 4.5 m above ground level.The sonic anemometer and the infrared CO 2 /H 2 O gas analyzer output were recorded at half-hour intervals by datalogger (CR5000, Campbell Scientific, Inc., Jackson, MS, USA), with high frequency of 10 Hz.Carbon and energy fluxes were calculated by using the half-hour covariance of vertical wind velocity and virtual temperature, water vapor density, and CO 2 density [15].We also measured some environmental factors that were used for gap filling calculations nearby the eddy tower.Photosynthetically active radiation (PAR) was observed by LI-190SB sensors (LI-COR Inc., Lincoln, NE, USA).Relative humidity and air temperature were observed by humidity/temperature probes (HMP45C, VAISALA Woburn, MA, USA).Precipitation was measured at the height of 1.5 m above the ground level by a tipping bucket rain-gauge (TE525MM, Campbell Scientific, Inc., Jackson, MS, USA).The copper-constantan thermocouples were installed at five depths (0.05, 0.10, 0.20, 0.50, and 1.0 m) below the ground to measure soil temperature.All data were saved in the digital datalogger (CR23X; Campbell Scientific, Inc., Jackson, MS, USA) [19].
For reducing the measurement-induced uncertainties, quality control was conducted for the observed data during 2009-2010.The CO 2 flux data observed during rainfall or snowfall or instrument failure (e.g., system maintenance, power outages, etc.) were eliminated as invalid data that has led to data gaps.To fill the CO 2 flux data gap, the gap filling methods, including mean diurnal variation (MDV), of previous periods and nonlinear regression were used in this study.For the MDV method, the means for different time interval (window size, usually 4-15 days) that were based on adjacent days were used to replace the missing flux data.Nonlinear regression method was interpolated by establishing the regression relationships between net ecosystem CO 2 exchanges (NEE) components and the associated environmental factors during the studied period [46].
Generally, the EC measured C fluxes represent NEE, which is the balance between ecosystem respiration (R e ) and GPP [47].Therefore, GPP can be obtained by the difference between R e and NEE.R e is summed by the daytime (R e,day ) and nighttime (R e,night ) ecosystem respiration.R e,night is derived from the nighttime NEE.Because R e,night is relevant to the soil temperature (T s ) [15,48], a temperature dependent model was derived from the observed R e,night with T s [19].Then, we estimated the R e,day rates and GPP by extrapolating the exponential relationship to the daytime periods.
R e,night = ae (bTs) where T s is soil temperature ( • C) at the depth of 0.05 m, and a and b are coefficients.During the data processing, the unit of CO 2 fluxes (NEE, GPP, R e , R e,day, R e,night ) is mg CO 2 m −2 s −1 .After that, we accumulated and converted all data to the whole day, and obtained the daily CO 2 flux data, when the unit of CO 2 fluxes (NEE, GPP, R e ) is g CO 2 m −2 day −1 or g C m −2 day −1 .
The daily GPP and climate data were aggregated to eight-day intervals in order to be consistent with MODIS eight-day composites.The aggregated eight-day GPP and climate data over 2009-2010 were utilized as parameter, input, and validated data in this study to support model simulation and validation.

Moderate Resolution Imaging Spectroradiometer Data and Vegetation Indices
The time series site-specific VIs for the targeted wetland EC tower were extracted from the MODIS datasets.We mainly used seven spectral bands to study vegetation and land surface: blue (459-479 nm), green (545-565 nm), red (620-670 nm), near infrared (NIR, 841-875 nm, 1230-1250 nm), and shortwave infrared (SWIR, 1628-1652 nm, 2105-2155 nm).In this study, based on the location (latitude and longitude) of the targeted EC tower, we downloaded the eight-day MOD09A1 data over 2009-2010 from Earth Observation and Modeling Facility (EOMF), the University of Oklahoma website (http://www.eomf.ou.edu).Then, we calculated EVI [37], LSWI [33], and NDVI [32] based on the reflectance values of the eight-day MODIS dataset from the blue, red, NIR, and SWIR spectral bands [33].Owing to some defects of the NDVI, the EVI was proposed by using a feedback-based approach that incorporates both background adjustment and atmospheric resistance concepts into the NDVI, resulting in a feedback-based, soil and atmosphere resistant vegetation index [37].Furthermore, VPM model used EVI and LSWI as important input data to simulate the dynamics of GPP at the targeted coastal wetland.
For time series of VIs calculated from surface reflectance, a simple approach was employed to fill those cloudy pixels [36,50].The gap-filling procedure that was employed for cloudy pixels of VIs was reported in the earlier studies [22,36,39].A three-point time series filter X(t − 1), X(t), and X(t + 1)) was employed to correct a cloudy pixel using values of non-cloudy pixels in the window.If X(t + 1) and X(t − 1) pixels were cloud-free, the average of X(t +1) and X(t − 1) was calculated and was used as X(t).If only one pixel (either X(t + 1) or X(t − 1)) was cloud-free, then we used that pixel to substitute X(t).If the three-point time series filter was not available, then we expanded it to a five-point time series filter (X(t − 2), X(t − 1), X(t), X(t + 1), and X(t + 2)) [36].The gap filling procedure was same as the above three-point time series filter.

Model Structure
As a RS-based vegetation photosynthesis model for simulating GPP [33], the VPM model is based on the conceptual partitioning of chlorophyll and non-photosynthetically active vegetation at the canopy.Air temperature (Ta), the EVI, LSWI, and PAR are the input data of the VPM model (Figure 2).
VPM estimates GPP over the photosynthetically active period of vegetation [33,36]: where parameter ε g is the light use efficiency (LUE) (µmol CO 2 /µmol photosynthetic photon flux density, PPFD), which is affected by leaf phenology, temperature, and water; parameter PAR is the photosynthetically active radiation (µmol PPFD), and parameter FPAR chl is the fraction of PAR absorbed by leaf chlorophyll in the canopy.
In the VPM model, FPAR chl is calculated as [10,23,33,36,39,40]: where EVI is the Enhanced Vegetation Index, which is calculated from the reflectance values of the eight-day MODIS dataset; α is the coefficient in the EVI-FPAR chl linear function, which is set to be 1.0 [33,36].ε g is calculated by the maximum light use efficiency (LUE) ε 0 (µmol CO 2 /µmol PPFD) and the scalars for the effects of temperature, water and leaf phenology on light use efficiency of vegetation.
T scalar is the scalar of temperature on photosynthesis, which is estimated by the Terrestrial Ecosystem Model (TEM) [51]; W scalar is the impact of water on vegetation photosynthesis with LSWI max [40]; and, P scalar is the scalar for the impact of leaf phenology on photosynthesis within the canopy [39].
For time series of VIs calculated from surface reflectance, a simple approach was employed to fill those cloudy pixels [36,50].The gap-filling procedure that was employed for cloudy pixels of VIs was reported in the earlier studies [22,36,39].A three-point time series filter X(t − 1), X(t), and X(t + 1)) was employed to correct a cloudy pixel using values of non-cloudy pixels in the window.If X(t + 1) and X(t − 1) pixels were cloud-free, the average of X(t +1) and X(t − 1) was calculated and was used as X(t).If only one pixel (either X(t + 1) or X(t − 1)) was cloud-free, then we used that pixel to substitute X(t).If the three-point time series filter was not available, then we expanded it to a five-point time series filter (X(t − 2), X(t − 1), X(t), X(t + 1), and X(t + 2)) [36].The gap filling procedure was same as the above three-point time series filter.

Model Structure
As a RS-based vegetation photosynthesis model for simulating GPP [33], the VPM model is based on the conceptual partitioning of chlorophyll and non-photosynthetically active vegetation at the canopy.Air temperature (Ta), the EVI, LSWI, and PAR are the input data of the VPM model (Figure 2).ε 0 , maximum light use efficiency; T scalar , W scalar , P scalar , the scalars for the effects of temperature, water and leaf phenology on light use efficiency of vegetation, respectively; FPAR chl , the fraction of PAR absorbed by leaf chlorophyll; PAR, the photosynthetically active radiation; LSWI max , the maximum Land Surface Water Index (LSWI); T max , T min , and T opt , the maximum, minimum, and optimal air temperature, respectively.

Model Parameterization
In order to normally run the VPM model for simulating GPP of coastal wetland using input data, such as climate observation and MODIS imagery, we conducted a series of model parameterization.
Different vegetation types have a different ε 0 value.It can be obtained from analysis of NEE and photosynthetic photon flux density of the EC tower site [52,53].In this study, we used the Michaelis-Menten function (Equation ( 13)) in order to estimate the ε 0 of local vegetation type by the nonlinear model between the observed NEE and incident PAR, according to the flux data at 30 min intervals.Here, we took 2016.67mg CO 2 mol −1 PAR as ε 0 during the study period: where GPP max (mg CO 2 m −2 s −1 ) is maximum gross primary production, R eco (mg CO 2 m −2 s −1 ) is ecosystem respiration, NEE (mg CO 2 m −2 s −1 ) is the net ecosystem CO 2 exchange, PAR is the photosynthetically active radiation (µmol PPFD).Based on the diurnal GPP flux data, we chose the maximum value of GPP (0.63 mg CO 2 m −2 s −1 ) during the study period as GPP max .Different vegetation types have different estimation values of T scalar , T min , T opt , and T max [33].The daily daytime mean temperature was calculated by using the daily maximum and minimum air temperature, rather than daily mean air temperature [33,54].In our study, if the air temperature of the study site falls below T min , T scalar is set to zero.Because T scalar is the scalar of temperature on plant photosynthesis activities, we consider the local meteorological conditions and plant growth (observed GPP) comprehensively.When considering the effect of low temperature on plant growth, we finally defined T min as 0 • C (cold damage to plants).Based on the observed maximum temperature during the two study years, we finally defined T max as 35 • C. We also defined site-specific T opt as 26 • C by exploring the correlation between the observed vegetation indices, GPP and air temperature over 2009-2010.
W scalar can be derived from the water-sensitive vegetation index LSWI and LSWI max .Site-specific LSWI max is estimated according to the time series MODIS data [26,33,36].As a water-sensitive vegetation index, LSWI max varies across years under different water condition.Using the observed eight-day MODIS reflectance data, we calculated the site-specific LSWI values over 2009-2010 and chose the maximum value during the growing season as LSWI max (0.27 and 0.26 for 2009 and 2010, respectively).

Model Evaluation
In order to test VPM's applicability for coastal wetland and to check the consistency between the modeled and the observed results, four statistical criteria were used for the model evaluation: (i) the coefficient of determination (R 2 ).Due to temporal autocorrelation within simulated GPP and EC-based GPP, generalized least square regression was conducted by using gls() function in R package nlme to adjust the residuals.Accordingly, the R 2  Cox and Snell was instead used as pseudo-R 2 by performing nagelkerke() function in package rcompanion; Equation (2), the root of mean square error (RMSE, Equation ( 14)); Equation (3), the relative mean deviation (RMD, Equation ( 15)); Equation (4), the relative error (RE) [55].The RMSE provides the prediction error of model by heavily weighting high errors.The RMD can smooth out the differences between the simulated and measured results by weighting all of the errors in the same way, which when close to 0 indicate the absence of bias for the model [15,56].We use SPSS version 20.0 (SPSS Inc., Chicago, IL, USA) and R package 3.2.1 [57] to conduct all of the statistical calculations.
where O i represent the EC estimated values and P i represent the model-simulated values.O is the mean of the EC estimated value and P is the mean of the model-predicted value.n is the number of observations.

Seasonal Dynamics of Hydrothermal Conditions, Vegetation Indices, and Gross Primary Production
The dynamics of the main meteorological variables are given in Figure 3.There was significant seasonal variation of daily air temperature (Ta) during 2009-2010, varying between −11.0 and 33.4 • C. The maximum and minimum values of Ta were found during July and January, respectively.But, no strong inter-annual difference among the different temperature regimes was found.The mean annual Ta were 13.9 and 13.3 • C in 2009 and 2010, respectively, both warmer than the long-term average (12.1 • C).The annual precipitation (rainfall and snowfall), measured at the tower in the experiment were 571.4 and 523.5 mm for 2009 and 2010, respectively, both higher than the long-term average (520.0mm).In 2009, there were 68 days with precipitation.In 2010, there were 66 days with precipitation (Figure 3).The results also show that the rainy season started in early June, and it ended in late August.
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 19 precipitation (Figure 3).The results also show that the rainy season started in early June, and it ended in late August.The EVI, NDVI, and LSWI curves can better demonstrate the development and the senescence of coastal wetland vegetation (Figure 4a-c).Both time series EVI and NDVI from 2009 to 2010 showed a strong seasonality for this coastal wetland ecosystem.They increased in spring and peaked in mid-August (21 August), and then declined and kept low in winter, corresponding well with the seasonal dynamics of GPPEC (Figure 4).However, seasonal variation in EVI within the growing seasons differed slightly with the NDVI in terms of phase and magnitude (Figure 4a,b).The values of NDVI were obviously higher than that of EVI from late March to early November in both years.The maximum values of EVI (0.21 and 0.32) were much lower than the peaks of NDVI (0.40 and 0.61) in 2009 and 2010, respectively.The LSWI also exhibited seasonal variation, but it did not vary with time as significantly, as did NDVI or EVI, especially during plant growing season (Figure 4c).As shown in Figure 4c, the peak LSWI values were 0.27 (28 July) and 0.26 (12 August) for 2009 and 2010, respectively, advancing about 8-24 days as compared with the peak values of NDVI and EVI.High EVI and LSWI in August indicated that the canopy was rapidly developed photosynthetically active, while the quick drop of both EVI and LSWI in late September indicated an early plant senescence (Figure 4).
The time series GPPEC also shows an obviously seasonal and inter-annual variability (Figure 4d).With the raise of Ta and PAR, the photosynthetic capacity of coastal wetland ecosystem gradually increased.GPPEC began to increase in mid-April (1 g C m −2 day −1 or higher), achieved the maximum values in early August (2009) or mid-July (2010), and it then decreased to below 1 g C m −2 day −1 in late October (2009) or early November (2010).The peak and annual values of GPPEC in 2009 were 9.05 g C m −2 and 1068.51g C m −2 , respectively; whereas, the respective values of GPPEC in 2010 were 9.88 g C m −2 and 1102.84 g C m −2 , respectively (Figure 4d).The C uptake period (CUP) of wetland plant, as defined by GPPEC > 1 g C m −2 day −1 , ranges from mid-April to early November, which can be partly explained by the variation of Ta and PAR (Figures 3 and 4d).Further, the growing season period for coastal wetland vegetation, as defined by Vis, agreed with the CUP defined by GPPEC over the two study years (Figure 4d).The EVI, NDVI, and LSWI curves can better demonstrate the development and the senescence of coastal wetland vegetation (Figure 4a-c).Both time series EVI and NDVI from 2009 to 2010 showed a strong seasonality for this coastal wetland ecosystem.They increased in spring and peaked in mid-August (21 August), and then declined and kept low in winter, corresponding well with the seasonal dynamics of GPP EC (Figure 4).However, seasonal variation in EVI within the growing seasons differed slightly with the NDVI in terms of phase and magnitude (Figure 4a,b).The values of NDVI were obviously higher than that of EVI from late March to early November in both years.The maximum values of EVI (0.21 and 0.32) were much lower than the peaks of NDVI (0.40 and 0.61) in 2009 and 2010, respectively.The LSWI also exhibited seasonal variation, but it did not vary with time as significantly, as did NDVI or EVI, especially during plant growing season (Figure 4c).As shown in Figure 4c, the peak LSWI values were 0.27 (28 July) and 0.26 (12 August) for 2009 and 2010, respectively, advancing about 8-24 days as compared with the peak values of NDVI and EVI.High EVI and LSWI in August indicated that the canopy was rapidly developed photosynthetically active, while the quick drop of both EVI and LSWI in late September indicated an early plant senescence (Figure 4).
The time series GPP EC also shows an obviously seasonal and inter-annual variability (Figure 4d).With the raise of Ta and PAR, the photosynthetic capacity of coastal wetland ecosystem gradually increased.GPP EC began to increase in mid-April (1 g C m −2 day −1 or higher), achieved the maximum values in early August (2009) or mid-July (2010), and it then decreased to below 1 g C m −2 day −1 in late October (2009) or early November (2010).The peak and annual values of GPP EC in 2009 were 9.05 g C m −2 and 1068.51g C m −2 , respectively; whereas, the respective values of GPP EC in 2010 were 9.88 g C m −2 and 1102.84 g C m −2 , respectively (Figure 4d).The C uptake period (CUP) of wetland plant, as defined by GPP EC > 1 g C m −2 day −1 , ranges from mid-April to early November, which can be partly explained by the variation of Ta and PAR (Figures 3 and 4d).Further, the growing season period for coastal wetland vegetation, as defined by Vis, agreed with the CUP defined by GPP EC over the two study years (Figure 4d).

Correlation between GPPEC, Vegetation Indices, and Air Temperature
The comparisons between GPPEC and VIs (EVI and NDVI) show that EVI has a better relationship with GPP in terms of the phase and amplitude (Figure 4).In order to further evaluate the biological significance of VIs in GPP predictions, the simple linear relationships between GPPEC and VIs were investigated (Figure 5).Correlation analyses show that the dynamics of GPPEC correlated well with that of the VIs (EVI, NDVI, and LSWI), both in 2009 and 2010 (Figure 5).However, EVI (R 2 = 0.65, p < 0.0001) has a better correlation with GPPEC than NDVI does (R 2 = 0.48, p < 0.0001).EVI can account for 65% of the variance in GPPEC, while NDVI only 48%.Above results indirectly prove that the improved vegetation index EVI has a higher predictive power to simulate the GPP of the coastal wetland ecosystem.The LSWI also has a strong relationship with GPPEC (R 2 = 0.41, p < 0.0001, Figure 5), showing that the land surface water content is important for the estimation of GPP during growing season.

Correlation between GPP EC , Vegetation Indices, and Air Temperature
The comparisons between GPP EC and VIs (EVI and NDVI) show that EVI has a better relationship with GPP in terms of the phase and amplitude (Figure 4).In order to further evaluate the biological significance of VIs in GPP predictions, the simple linear relationships between GPP EC and VIs were investigated (Figure 5).Correlation analyses show that the dynamics of GPP EC correlated well with that of the VIs (EVI, NDVI, and LSWI), both in 2009 and 2010 (Figure 5).However, EVI (R 2 = 0.65, p < 0.0001) has a better correlation with GPP EC than NDVI does (R 2 = 0.48, p < 0.0001).EVI can account for 65% of the variance in GPP EC , while NDVI only 48%.Above results indirectly prove that the improved vegetation index EVI has a higher predictive power to simulate the GPP of the coastal wetland ecosystem.The LSWI also has a strong relationship with GPP EC (R 2 = 0.41, p < 0.0001, Figure 5), showing that the land surface water content is important for the estimation of GPP during growing season.
We also analyzed the relationships between VIs (EVI, NDVI, and LSWI), GPP EC , and air temperature (Ta).The regression analysis shows that Ta was positively correlated with NDVI (R 2 = 0.55, p < 0.0001), EVI (R 2 = 0.62, p < 0.0001), and GPP EC (R 2 = 0.78, p < 0.0001), respectively (Figure 6a,b,d).NDVI, EVI, and GPP EC gradually increased as the temperature rose, and the temporal variation in NDVI, EVI, and GPP EC fluxes could be largely explained by air temperature.LSWI had a significant U-shaped relationship (R 2 = 0.28, p < 0.01) with air temperature.As temperature rose, LSWI decreased first, reached its bottom, and then increased.The demarcation point is about 10 • C (Figure 6c).has a better correlation with GPPEC than NDVI does (R 2 = 0.48, p < 0.0001).EVI can account for 65% of the variance in GPPEC, while NDVI only 48%.Above results indirectly prove that the improved vegetation index EVI has a higher predictive power to simulate the GPP of the coastal wetland ecosystem.The LSWI also has a strong relationship with GPPEC (R 2 = 0.41, p < 0.0001, Figure 5), showing that the land surface water content is important for the estimation of GPP during growing season.

Simulation and Evaluation of Vegetation Photosynthesis Model
We run the VPM model based on the observation of Ta, PAR, and MODIS VIs (EVI and LSWI), and examined the seasonal variation of the modeled GPP (GPPVPM) for the tested coastal wetland We also analyzed the relationships between VIs (EVI, NDVI, and LSWI), GPPEC, and air temperature (Ta).The regression analysis shows that Ta was positively correlated with NDVI (R 2 = 0.55, p < 0.0001), EVI (R 2 = 0.62, p < 0.0001), and GPPEC (R 2 = 0.78, p < 0.0001), respectively (Figure 6a,b,d).NDVI, EVI, and GPPEC gradually increased as the temperature rose, and the temporal variation in NDVI, EVI, and GPPEC fluxes could be largely explained by air temperature.LSWI had a significant U-shaped relationship (R 2 = 0.28, p < 0.01) with air temperature.As temperature rose, LSWI decreased first, reached its bottom, and then increased.The demarcation point is about 10 °C (Figure 6c).

Simulation and Evaluation of Vegetation Photosynthesis Model
We run the VPM model based on the observation of Ta, PAR, and MODIS VIs (EVI and LSWI), and examined the seasonal variation of the modeled GPP (GPPVPM) for the tested coastal wetland during 2009-2010 (Figure 7).GPPVPM was close to zero from November to early March, then began to increase rapidly in mid-March, and achieved a peak in early-August and mid-July in 2009 and 2010, respectively.GPPVPM declined gradually after reaching peak it and dropped below 1 g C m −2 day −1 by early-November.Comparisons between GPPVPM and GPPEC fluxes at eight-day time intervals show that the simulated eight-day GPPVPM fluxes were in good agreement with GPPEC fluxes regarding to

Simulation and Evaluation of Vegetation Photosynthesis Model
We run the VPM model based on the observation of Ta, PAR, and MODIS VIs (EVI and LSWI), and examined the seasonal variation of the modeled GPP (GPP VPM ) for the tested coastal wetland during 2009-2010 (Figure 7).GPP VPM was close to zero from November to early March, then began to increase rapidly in mid-March, and achieved a peak in early-August and mid-July in 2009 and 2010, respectively.GPP VPM declined gradually after reaching peak it and dropped below 1 g C m −2 day −1 by early-November.Comparisons between GPP VPM and GPP EC fluxes at eight-day time intervals show that the simulated eight-day GPP VPM fluxes were in good agreement with GPP EC fluxes regarding to the patterns and magnitudes.VPM fairly reappeared the U-shaped growth pattern, which driveled by the variation in Ta and PAR, and accurately simulated the seasonal dynamics of GPP EC fluxes for the simulated two years.However, eight-day GPP VPM peak values were slightly higher than GPP EC in both 2009 and 2010 (Figure 7). the patterns and magnitudes.VPM fairly reappeared the U-shaped growth pattern, which driveled by the variation in Ta and PAR, and accurately simulated the seasonal dynamics of GPPEC fluxes for the simulated two years.However, eight-day GPPVPM peak values were slightly higher than GPPEC in both 2009 and 2010 (Figure 7).To further evaluate the accuracy of model simulation, a linear regression analysis between GPPVPM and GPPEC was conducted (Figure 8).For both years, the linear regression of GPPVPM against GPPEC shows a good correlation (R 2 ≥ 0.72, p < 0.0001; Figure 8; Table 1).The RMD and RMSE values were −1.01%and 25.29% during 2009-2010, respectively.The simulated annual GPPVPM were 1057.64 and 1091.76 g C m −2 yr −1 for 2009 and 2010, respectively, which is slightly lower than the estimated annual GPPEC (1068.51 and 1102.84 g C m −2 yr −1 for 2009 and 2010, respectively), both with relative error (RE) values of −1% (Table 1).Moreover, cumulative GPP fluxes simulated over the two years were in good agreement with the estimated values, with −1% RE (Table 1).To further evaluate the accuracy of model simulation, a linear regression analysis between GPP VPM and GPP EC was conducted (Figure 8).For both years, the linear regression of GPP VPM against GPP EC shows a good correlation (R 2 ≥ 0.72, p < 0.0001; Figure 8; Table 1).The RMD and RMSE values were −1.01%and 25.29% during 2009-2010, respectively.The simulated annual GPP VPM were 1057.64 and 1091.76 g C m −2 yr −1 for 2009 and 2010, respectively, which is slightly lower than the estimated annual GPP EC (1068.51 and 1102.84 g C m −2 yr −1 for 2009 and 2010, respectively), both with relative error (RE) values of −1% (Table 1).Moreover, cumulative GPP fluxes simulated over the two years were in good agreement with the estimated values, with −1% RE (Table 1).
Table 1.The evaluation of model simulation results in the typical coastal wetland.Four statistical criteria for the validation tests: pseudo-R 2 , the coefficient of determination; RMSE, the root of mean square error; RMD, the relative mean deviation; and, RE, the relative error.1).Moreover, cumulative GPP fluxes simulated over the two years were in good agreement with the estimated values, with −1% RE (Table 1).

Biophysical Performance of Vegetation Indices in Typical Coastal Wetland
Monitoring vegetation change over time at larger scales using vegetation indices has greatly enhanced our understanding of the changing environment [58,59].Different vegetation indices vary in their biophysical expression and they can offer distinct information for the leaf phenological cycle, which in turn, affect plant photosynthesis and regulate ecosystem carbon balance [10,22,36].Some researchers have reported that time series of NDVI, EVI (greenness-related vegetation indices), and LSWI (water-related vegetation index) can be used to delineate the green-up and the senescence phase at the canopy level [36,60].However, which is the best vegetation index for the model to more accurately predict the plant production in coastal wetlands remains unknown.Therefore, it is important to evaluate the biophysical and phenological expression of different vegetation indices in coastal wetland ecosystems to accurately assess vegetation phenology and predict plant production [61,62].
In this study, we examined the relationship between NDVI, EVI, LSWI, and estimated GPP EC , Ta, and evaluated the biophysical expression of vegetation indices in a typical coastal wetland.Both EVI and NDVI have significant correlations with GPP EC , but EVI has a stronger correlation with GPP EC , and hence has more biological and phenological significance in GPP predictions than the NDVI in this coastal wetland (Figures 5 and 6).EVI can minimize the influences of soil condition and atmosphere in reflectance data by incorporating soil adjustment and atmospheric resistance factors [58].Time series of LSWI also provide valuable information (e.g., water status) for simulating plant growth and wetland carbon exchange, even though it has a lower relationship between Ta and GPP EC .It has a great potential in obtaining the water content information of leaf and canopy, and can more accurately identify the vegetation phenology by combining both NIR and SWIR bands [22,63,64].The results in this study are consistent with earlier studies [10,15,22,23,33,40,65,66], in which all indirectly prove the hypothesis of chlorophyll-FPAR chl -EVI and leaf 1 water-LSWI, as applied in the VPM model.However, it is necessary to further develop the LSWI-based phenology algorithm in the near future [40].

Model Comparison and Error Source Analysis
Recently, using the LUE models to simulate GPP plays an important role in the global carbon budget research, but our understanding of relative advantages and disadvantages of different models is still limited [67].The VPM model simulated the trends of the coastal wetland phenology well and precisely predicted GPP (R 2 = 0.73) in this study.The MODIS-Photosynthesis (PSN) model also can simulate the global GPP and generate the Terra/MODIS GPP project (MOD17A2) [40,68].Here, eight-day MODIS GPP product (GPP MOD17A2 ) data during 2009-2010 were downloaded from the Oak Ridge National Laboratory's Distributed Active Archive Center (DAAC) website (http://daac.ornl.gov/MODIS)[69,70].Then, we compared the differences between GPP EC , GPP VPM , and GPP MOD17A2 in our studied coastal wetland.A Taylor diagram was also used to display the quality of model predictions against the observations [71].We found that both GPP EC and GPP VPM were substantially higher than GPP MOD17A2 in terms of the magnitudes.The underestimated range of GPP MOD17A2 is probably 40-80% (see Figure 9).Some previous studies also reported that the MODIS-PSN model underestimated GPP in different ecosystems [15,40,[72][73][74].However, Jin et al., found that GPP MOD17A2 was significantly lower than GPP EC during the prophase of the growing season, but higher than GPP EC during the late growing season at broadleaf deciduous woodland [39].
The discrepancies between GPPMOD17A2 and GPPEC may result from two underestimated sources.First, GPPMOD17A2 uses global meteorological reanalysis dataset, including average and minimum Ta, incident PAR, and specific humidity.The input weather data of MODIS-PSN model may not match the local air temperature data.The second source is the parameter ofε0 in MODIS-PSN [66,75].ε0 is a very important parameter in the PEM models to estimate GPP and differs significantly among vegetation types [76].However, the ε0 and some other biome-specific physiological parameters are not differentiated from the different performance of a specific biome [77].The EC technique provides an appropriate method to accurately estimate the parameter ε0 by using a large number of measured data [23,30].In this study, we estimated the ε0 value by fitting the nonlinear model between observed NEE and PAR in 2009 and 2010.Then, we used the site-specific climate data, the individual ε0 estimated from CO2 flux data, and vegetation indices (EVI, LSWI) as input data and the parameters of VPM to simulate the GPP dynamics of coastal wetland.Therefore, future works should pay more attention to the local climate data and biome-specific physiological parameters in using GPPMOD17A2 for regional analysis.

Sensitivity and Uncertainty for Vegetation Photosynthesis Model Simulations
In general, there are always some errors in model simulations.The analysis of sources and uncertainties for the errors are very important for future improvement and the development of these models [78,79].In this study, the VPM simulation results show that the VPM model can correctly simulate the seasonal variations of GPP.The modeled results can explain about 85% of GPPEC changes over the two study years.However, there still exist some differences between seasonal GPPVPM and GPPEC in this study, such as higher GPPVPM from late July to early August in 2009 and 2010, and lower The discrepancies between GPP MOD17A2 and GPP EC may result from two underestimated sources.First, GPP MOD17A2 uses global meteorological reanalysis dataset, including average and minimum Ta, incident PAR, and specific humidity.The input weather data of MODIS-PSN model may not match the local air temperature data.The second source is the parameter ofε 0 in MODIS-PSN [66,75].ε 0 is a very important parameter in the PEM models to estimate GPP and differs significantly among vegetation types [76].However, the ε 0 and some other biome-specific physiological parameters are not differentiated from the different performance of a specific biome [77].The EC technique provides an appropriate method to accurately estimate the parameter ε 0 by using a large number of measured data [23,30].In this study, we estimated the ε 0 value by fitting the nonlinear model between observed NEE and PAR in 2009 and 2010.Then, we used the site-specific climate data, the individual ε 0 estimated from CO 2 flux data, and vegetation indices (EVI, LSWI) as input data and the parameters of VPM to simulate the GPP dynamics of coastal wetland.Therefore, future works should pay more attention to the local climate data and biome-specific physiological parameters in using GPP MOD17A2 for regional analysis.

Sensitivity and Uncertainty for Vegetation Photosynthesis Model Simulations
In general, there are always some errors in model simulations.The analysis of sources and uncertainties for the errors are very important for future improvement and the development of these models [78,79].In this study, the VPM simulation results show that the VPM model can correctly simulate the seasonal variations of GPP.The modeled results can explain about 85% of GPP EC changes over the two study years.However, there still exist some differences between seasonal GPP VPM and GPP EC in this study, such as higher GPP VPM from late July to early August in 2009 and 2010, and lower GPP VPM in late May in 2010 (Figure 7).The discrepancies between GPP EC and GPP VPM can be attributed, in part, to the estimation error of tower-based GPP EC [36,80].GPP EC is calculated as the difference between ecosystem respiration and NEE [15,36].The estimation error of daytime ecosystem respiration, which is the difference between ecosystem respiration and observed nighttime ecosystem respiration, may lead to the error of estimation (either overestimation or underestimation) of GPP EC [15,36,39,40].
Another error source may be attributed to the modeling of GPP VPM from VPM [33].To further clarify the sources and uncertainty of the GPP VPM errors, we conducted a sensitivity analysis for some parameters and variables of the model (Figure 10).The alternative sensitivity scenarios of input parameter were set by increasing Ta, ε 0 , PAR, EVI, and LSWI by 10%, respectively.When compared with the base scenario, the sensitivity of ε 0 and PAR are relatively larger than that of Ta in this coastal wetland (Figure 10a,d).The potential conversion efficiency of absorbed photosynthetically active radiation is largely determined by ε 0 , which varies with vegetation types and study regions [23,39], and it may even affect the accuracy of GPP simulation [40,81].PAR is an important variable in the GPP simulation, lower/higher PAR values may result in the under/over-estimation of GPP for the assumed relationship between GPP and PAR in VPM model [15,82].Therefore, the accurate estimation of ε 0 and development of PAR observation would substantially improve the simulated accuracy of the VPM model and other models [36].The source of uncertainty also come from the MODIS time series vegetation indices, such as the effects of cloud, cloud shadow, aerosols, and other atmospheric condition [40].The EVI is more sensitive than LSWI to GPP simulation during 2009-2010 in this typical coastal wetland (Figure 10b,d).Future work should focus on how to avoid and mitigate the effects of angular geometry on surface albedo and to reconstruct the time series data of vegetation indices in order to improve the accuracy of vegetation indices, especially EVI.
The simulation uncertainties of VPM also come from two intermediate variables that are related to water (W scalar ) and temperature (T scalar ) [39].We use the two intermediate variables to set three sensitivity scenarios, including without T scalar (ε g = W scalar × P scalar × ε 0 ), without W scalar (ε g = T scalar × P scalar × ε 0 ), and without T scalar & W scalar (ε g = P scalar × ε 0 ).When there is no T scalar in the VPM, the model overestimated GPP by about 10%, 39%, and 25% in 2009, 2010, and 2009-2010, respectively.In the absence of W scalar in VPM, the model overestimated GPP by about 18%, 36%, and 27% in 2009, 2010, and 2009-2010, respectively.When there is no both T scalar and W scalar in VPM, the model overestimated GPP by about 31%, 57%, and 44% in 2009, 2010, and 2009-2010, respectively (Figure 10c,d).Above results indicated that both T scalar and W scalar have profound effects on the predicted accuracy of GPP.It will be essential to integrate both temperature and water to downscale maximum light use efficiency to better represent the impacts of temperature and water on plant photosynthesis and GPP [36,39,40,66].

Conclusions
In this study, we incorporated climate and remote sensing data into a satellite-based VPM model to investigate the biophysical performance of different vegetation indices and to evaluate the applicability and accuracy of the VPM model for simulating GPP of coastal wetland in China.The results indicated that air temperature can explain most of the seasonal variation in vegetation indices and GPP EC fluxes.Time series of greenness-related EVI and water-related LSWI can better delineate the green-up and senescence phases of coastal wetland vegetation, which are in accordance with the CUP, as defined by GPP EC .Further, both EVI and water-sensitive LSWI have higher predictive power for simulating the GPP in this coastal wetland.The EVI-based VPM model is able to simulate the U-shaped growth pattern of the coastal plant and accurately simulate the trajectories of the seasonal dynamics of GPP EC fluxes, explaining about 85% of GPP EC changes.Further development and improvement of VPM for nearby coastal wetlands and the spatial-temporal scale expansion are necessary for regional and continental simulations of GPP in the future.

Figure 1 .
Figure 1.Study site location at the Yellow River Delta.The remote sensing data is derived from Google Earth.The eddy covariance system, including a flux tower and a weather station, was deployed at the coastal zone of Yellow River Delta Nature Reserve in 2008, with observation of CO2 and water fluxes, air temperature, precipitation, photosynthetically active radiation (PAR), etc.

Figure 1 .
Figure 1.Study site location at the Yellow River Delta.The remote sensing data is derived from Google Earth.The eddy covariance system, including a flux tower and a weather station, was deployed at the coastal zone of Yellow River Delta Nature Reserve in 2008, with observation of CO 2 and water fluxes, air temperature, precipitation, photosynthetically active radiation (PAR), etc.

Figure 2 .Figure 2 .
Figure 2. Structure of the vegetation photosynthesis model (VPM) model.CO2 flux and climate data (air temperature, PAR, precipitation, etc.) come from the eddy flux tower.Vegetation indices (VIs) were extracted from the Moderate Resolution Imaging Spectroradiometer (MODIS) datasets.ε0, Figure 2. Structure of the vegetation photosynthesis model (VPM) model.CO 2 flux and climate data (air temperature, PAR, precipitation, etc.) come from the eddy flux tower.Vegetation indices (VIs) were extracted from the Moderate Resolution Imaging Spectroradiometer (MODIS) datasets.ε0 , maximum light use efficiency; T scalar , W scalar , P scalar , the scalars for the effects of temperature, water and leaf phenology on light use efficiency of vegetation, respectively; FPAR chl , the fraction of PAR absorbed by leaf chlorophyll; PAR, the photosynthetically active radiation; LSWI max , the maximum Land Surface Water Index (LSWI); T max , T min , and T opt , the maximum, minimum, and optimal air temperature, respectively.

Figure 3 .
Figure 3. Seasonal variation in daily air temperature, precipitation, and aggregated eight-day PAR from eddy covariance (EC) tower during 2009-2010.

Figure 3 .
Figure 3. Seasonal variation in daily air temperature, precipitation, and aggregated eight-day PAR from eddy covariance (EC) tower during 2009-2010.

Figure 7 .
Figure 7. Comparisons of the seasonal dynamics between the estimated gross primary production from EC tower (GPPEC) and simulated gross primary production by VPM (GPPVPM) during 2009-2010.Spring: March to May; Summer: June to August; Autumn: September to November; and, Winter: December to February.

Figure 9 .
Figure 9. Comparisons between the estimated gross primary production from the flux tower (GPPEC) and simulated gross primary production by VPM (GPPVPM) and MODIS global GPP data product (GPPMOD17A2) in the typical coastal wetland during 2009-2010.Taylor diagram provides a visual comparison between two model results, or, most commonly, between model results and observations [71].

Figure 9 .
Figure 9. Comparisons between the estimated gross primary production from the flux tower (GPP EC ) and simulated gross primary production by VPM (GPP VPM ) and MODIS global GPP data product (GPP MOD17A2 ) in the typical coastal wetland during 2009-2010.Taylor diagram provides a visual comparison between two model results, or, most commonly, between model results and observations [71].