Estimation and Analysis of Spatiotemporal Dynamics of the Net Primary Productivity Integrating Efficiency Model with Process Model in Karst Area

Estimates of regional net primary productivity (NPP) are useful in modeling regional and global carbon cycles, especially in karst areas. This work developed a new method to study NPP characteristics and changes in Chongqing, a typical karst area. To estimate NPP accurately, the model which integrated an ecosystem process model (CEVSA) with a light use efficiency model (GLOPEM) called GLOPEM-CEVSA was applied. The fraction of photosynthetically active radiation (fPAR) was derived from remote sensing data inversion based on moderate resolution imaging spectroradiometer atmospheric and land products. Validation analyses showed that the PAR and NPP values, which were simulated by the model, matched the observed data well. The values of other relevant NPP models, as well as the MOD17A3 NPP products (NPP MOD17), were compared. In terms of spatial distribution, NPP decreased from northeast to southwest in the Chongqing region. The annual average NPP in the study area was approximately 534 gC/m2a (Std. = 175.53) from 2001 to 2011, with obvious seasonal variation characteristics. The NPP from April to October accounted for 80.1% of the annual NPP, while that from June to August accounted for 43.2%. NPP changed with the fraction of absorbed PAR, and NPP was also significantly correlated to precipitation and temperature at monthly temporal scales, and showed stronger sensitivity to interannual variation in temperature.


Introduction
Net primary productivity (NPP) is the net flux of carbon from the atmosphere into green plants per unit time.It is an essential flux of the net ecosystem exchange of CO 2 between the atmosphere and terrestrial ecosystems.It is a key component of the carbon cycle [1] and an important indicator of the ecosystem [2].Regional estimates of NPP are useful in modeling the regional and global carbon cycle.Quantitative estimates of NPP at regional to global scales are significant for understanding changes in ecosystem structure and function, predicting terrestrial carbon cycle trends [3][4][5], determining sustainable use of natural resources, and making policy decisions [6].
Three main methods have been used to estimate NPP for regional areas: (1) the establishment of a model that combined climatic and environmental factors, was established to evaluate the regional distribution of carbon budgets [7,8]; (2) the application of process-based, biogeochemical models that scale-up and extrapolate site-specific measurements [9,10]; and (3) models that use remote sensing data, numerous efforts have been made to develop various NPP models and to estimate global NPP.The accuracy of satellite-derived NPP estimates have been recently highlighted as a matter of concern because of poor satellite data quality and lack of accurate ecosystem representation in satellite-based models [11][12][13].
Numerous remote sensing-based models for estimating NPP, such as Carnegie Ames Stanford Approach (CASA) model [14], Simulation model of Carbon Cycle in Land Ecosystems (MOD-Sim-Cycle) model [15], and GLOBAL Production Efficiency Model (GLOPEM) model [16], have been developed in recent decades to study the dynamics of vegetation productivity.Examples of ecological process models that simulate basic ecological physiology and the process of plant growth include Boreal Ecosystem Productivity Simulator (BEPS) model, Terrestrial Ecosystem Model (TEM) model, CENTURY model, and BioGeochemical Cycles (Biome-BGC) model [17].The production efficiency model based on satellite remote sensing encounters difficulty in explaining the response of ecosystem internal mechanism to environmental change.Affected by the method and the data quality, it will produce deviation in analysis of large scale ecosystem changes.Thus, an integrated strategy was developed that combined the remote sensing-based GLOBAL Production Efficiency Model (GLOPEM) and an ecosystem process-based model, the Carbon Exchange between Vegetation, Soil and Atmosphere model (CEVSA), enhancing the possibility and feasibility of NPP model estimation.
The GLOPEM-CEVSA model [18] makes full use of the advantages of remote sensing and considers the spatial heterogeneity of ecosystems.The model is suitable for estimating the spatial-temporal distribution of NPP for terrestrial ecosystems.As the remote sensing data have the timeliness and space, for the model, the time and space scale can be changed, this paper selected 2001-2011, a total of 11 years of remote sensing products.However, few available models integrate moderate-resolution imaging spectroradiometer (MODIS) atmosphere and land products to estimate fraction of photosynthetically active radiation (fPAR).A main component of the GLOPEM-CEVSA modeling framework is the satellite-derived FPAR absorbed by vegetation, which can be derived from satellite normalized difference vegetation index (NDVI) data that come from MODIS products, with a 1000 m resolution.
The LUE model is an index of the photosynthesis.Research of LUE is beneficial for the energy cycle of the biosphere and of great significance for energy use and transformation, and the mechanism of vegetation physiological and ecological processes.Considering the atmosphere, radiation, moisture, carbon dioxide and other stress factors as a model-driven parameter, the model is used to simulate the autotrophic respiration of plants and calculate the dynamic NPP.Furthermore, soil moisture data were obtained by inversion of precipitation data and evapotranspiration models.In the southwest region of China, especially the Chongqing region, soil is shallow, and the proportion of vegetation carbon to carbon reservoir is much higher than that of other types.Understanding carbon sources and sinks is crucial to provide scientific reference for the management of rocky desertification in karst forest."karst" refers to a type of terrain, usually formed on carbonate rock where groundwater has solutionally enlarged openings to form a subsurface drainage system.Most of the studies on carbon cycle in karst are focused on the micro scale of ecosystem, with few studies on regional scale.NPP provides the energetic and material basis for all heterotrophic life.Thus, there is a need to study problems from the perspective of carbon cycles and energy flows.
The main aim of this article is to estimate and analyze NPP in the Chongqing area by using GLOPEM-CEVSA model with MODIS land parameters.The objectives are as follows: (1) to optimize of GLOPEM-CEVSA model parameters, invert fPAR using remote sensing data, and adapt and validate a remote sensing-driven parametric model to calculate photosynthetically active radiation (PAR); (2) to investigate the seasonal variation and spatial distribution of NPP from 2001 to 2011 in the Chongqing area; and (3) to analyze how the climate and vegetation factors influence NPP.

Study Area
Chongqing is located in a typical karst rocky desertification region of southwestern China (105   north-south dimension of 450 km, and is 470 km from east to west.It is in the upper reaches of the Yangtze River valley.Chongqing (Figure 1) lies in the transition zone of the Qinghai-Tibet Plateau and the Yangtze River valley in the middle and lower reaches, and is administered by 38 counties (autonomous counties).Three major tectonic units are located in the area, namely, the Daba Mountain fault fold belt, the eastern Sichuan fold belt and the Sichuan-Hubei-Hunan-Guizhou uplift belt.The main landforms are mountains and hills.Most of the region is located in the karst restoration areas.The Yangtze River from west to east stretches across an area of 665 km.Chongqing has a humid subtropical monsoon climate with an annual average temperature of 16-18 • C, an average winter temperature of 6-8 • C, and abundant rainfall.The weather in the second half of the year is cloudy with short days and the average annual precipitation is 1300 mm.Chongqing is mostly karst area, with only a small part of the northwest being non-karst area.Karst ecosystems are very fragile.Analysis of the productivity of vegetation, assessment of the carbon/sink in the region quantitatively, and understanding the change pattern and trend of the vegetation ecological effect are of great significance to carry out ecological restoration.
south dimension of 450 km, and is 470 km from east to west.It is in the upper reaches of the Yangtze River valley.Chongqing (Figure 1) lies in the transition zone of the Qinghai-Tibet Plateau and the Yangtze River valley in the middle and lower reaches, and is administered by 38 counties (autonomous counties).Three major tectonic units are located in the area, namely, the Daba Mountain fault fold belt, the eastern Sichuan fold belt and the Sichuan-Hubei-Hunan-Guizhou uplift belt.The main landforms are mountains and hills.Most of the region is located in the karst restoration areas.The Yangtze River from west to east stretches across an area of 665 km.Chongqing has a humid subtropical monsoon climate with an annual average temperature of 16-18 °C, an average winter temperature of 6-8 °C, and abundant rainfall.The weather in the second half of the year is cloudy with short days and the average annual precipitation is 1300 mm.Chongqing is mostly karst area, with only a small part of the northwest being non-karst area.Karst ecosystems are very fragile.Analysis of the productivity of vegetation, assessment of the carbon/sink in the region quantitatively, and understanding the change pattern and trend of the vegetation ecological effect are of great significance to carry out ecological restoration.

Remote Sensing Data
This study used MODIS NDVI data derived from MOD13A3_V005 (Grade 3 NDVI Vegetation Index Products) with a spatial resolution of 1 km and a temporal resolution of 1 month, covering the period of January 2001 to December 2011.The Chongqing area includes h27v05 and h27v06 strips.The data are mainly used to invert fPAR.All these three types of data came from the NASA-MODIS website (http://ladsweb.nascom.nasa.gov)and were openly accessible.The MODIS/NDVI 1 km products were extracted from the MODIS/NDVI data [19] using MODIS Reprojection Tools software.The 264 scene images were synthesized and mosaicked, and the image format and projection were converted.The MODIS stored SIN coordinates were converted to WGS84.The map projection was converted to Lambert conformal conic projection, converting the format of MODIS remote sensing images from HDF scientific data format to TiFF.
The latest correction method and the correction coefficient for radiation correction are based on remote sensing data.The measured GPS control points on the satellite image for geometric correction are used to ensure that the precision of geometric correction is better than ponepixel.Atmospheric radiation correction and reflectance transformation are conducted using linear transformation.

Remote Sensing Data
This study used MODIS NDVI data derived from MOD13A3_V005 (Grade 3 NDVI Vegetation Index Products) with a spatial resolution of 1 km and a temporal resolution of 1 month, covering the period of January 2001 to December 2011.The Chongqing area includes h27v05 and h27v06 strips.The data are mainly used to invert fPAR.All these three types of data came from the NASA-MODIS website (http://ladsweb.nascom.nasa.gov)and were openly accessible.The MODIS/NDVI 1 km products were extracted from the MODIS/NDVI data [19] using MODIS Reprojection Tools software.The 264 scene images were synthesized and mosaicked, and the image format and projection were converted.The MODIS stored SIN coordinates were converted to WGS84.The map projection was converted to Lambert conformal conic projection, converting the format of MODIS remote sensing images from HDF scientific data format to TiFF.
The latest correction method and the correction coefficient for radiation correction are based on remote sensing data.The measured GPS control points on the satellite image for geometric correction are used to ensure that the precision of geometric correction is better than ponepixel.Atmospheric radiation correction and reflectance transformation are conducted using linear transformation.

Meteorological Data
The meteorological data, derived from 34 meteorological stations in Chongqing from 2001 to 2011, were provided by the climate center of Chongqing Meteorological Bureau.The data include monthly average temperature, monthly total precipitation, wind speed, sunshine hours, and relative humidity.A meteorological raster map with a spatial resolution of 1 km and a temporal resolution of 1 month was obtained by interpolating meteorological data from the ANUSPLIN [20], an interpolation software developed by the Australian National University using a smooth thin-plate spline method to enter.

Model Parameter Data
The model-assisted parameter data (Table 1) mainly include the Chongqing area administrative division map, latitude and longitude, altitude, land cover data, and digital elevation model (DEM), among others.Land cover data are derived from a LANDSAT-TM image with the aid of MCD12Q1 [21].MCD12Q1 is from the Earth Observation System/Medium resolution Imaging Spectroradiometer (EOS/MODIS) of the National Aeronautics and Space Administration (NASA), with a spatial resolution of 500 m × 500 m from 2001 to 2009.The vegetation biomass data are the statistical of the area of the dominant tree species provided by the National Forest Resources Statistics (1989)(1990)(1991)(1992)(1993).According to the relationship between the biomass (B) and the storage volume (V), the total biomass of each dominant tree species and the average biomass of each forest vegetation type were calculated.The results of the study on the biomass of forest, farmland, grassland and desert in China obtained by Fang et al. (1996) were referenced.In addition, the values given by the CEVSA model, and the biomass allocation parameters of each stand refer to the values given by the Biome-BGC model and GLOPEM-CEVSA model that were improved by Wang Junbang.The DEM data with a spatial resolution of 1 km were obtained by spatial resampling using 90 m resolution STRM data provided by the Digital Image Network of the Computer Network Information Center of the Chinese Academy of Sciences.After removing outliers, the projection and resampling processes are redefined to be consistent with the rest of the study data presented in this paper.First, the NPP is calculated from GPP minus autotrophic respiration (Ra) based on light energy efficiency model using fPAR, PAR and meteorological data.The relevant parameters of Ra can be obtained from the Biome-BGC ecological process model.The reliability of the results is verified by comparing the simulation results of GLOPEM-CEVSA with those of the other relevant NPP models in Chongqing.The model is validated in two ways.The first uses forest inventory data, in which the NPP was verified by the biomass conversion factor method.The converted value corresponds to the result in time and space.This precision testing is the most convincing method.The second is a cross-validation among the relevant models, especially the MODIS NPP product MOD17.

GLOPEM Model
GlOPEM model [22][23][24] utilizes the production efficiency concept, in which the canopy absorbed photosynthetically active radiation (APAR) is used with a conversion "efficiency".The model is simulated by an empirical relationship between vegetation biomass and air temperature, which does not distinguish plant organs, or include biomass allocation simulation.
where GPP is Gross Primary Productivity (GPP) over given time intervalt, fPAR is fraction of absorbed photosynthetically active radiation, ∈ m is maximum light use efficiency and f(T) and f(W) are the constraints resulting from temperature to water stress.Ra is the autotrophic respiration, w is the above ground biomass, Tc is climatological average air temperature, T is air temperature, and ρ min is minimum value of the reflectivity in visible light channel.The NPP is further estimated by adjusting autotrophic respiration from GPP.

CEVSA Model
The CEVSA model is the vegetation submodel estimating vegetation distribution and calculating net primary production (NPP) [25], carbon allocation and litter production, which simulates synthesis of the biophysical and biogeochemical processes between vegetation, soil and the atmosphere.
The model consists of three sub-modules, biophysical module, plant physiological growth module, and soil carbon reservoir module.Biophysical module is primarily used to calculate the hydrothermal exchanges between atmosphere and land, and the dynamics of soil moisture.Plant physiological growth module is mainly used to calculate the accumulation and distribution of carbon and nitrogen between various organs of plants, as well as the plant photosynthesis and respiration.Soil carbon reservoir module is used to estimate the decomposition and transformation of soil organic matter.
where A b and A d are the assimilation rate determined by enzyme system and stomatal conductance; g s is stomatal conductance to water vapor; W c is the carboxylation rate; W j is the carboxylation rate depended on the rate of electron transport; W p is the carboxylation rate; P 0 and P c are the internal partial pressures of O 2 and CO 2 , respectively; τ is the specificity factor for CO 2 ; R d is the rate of respiration in light; Pa is the partial pressure of CO 2 in the leaf surface; R h is relative humidity of the air surrounding the leaf; K g (w s ) is a hyperbolic response function which describes the influence of soil water content w s on stomatal conductance g s ; c p is the heat of air; g n is the stomatal conductance; g a is the boundary layer conductance; R n is the net radiation; γ is the psychrometric constant; σ is the heat of vaporization; ρ is the mean air density; D is vapor pressure deficit; s is the slope of saturation vapour pressure; w s is the soil water content; and Rain, ET and Runoff are the rainfall reaching the ground, evapotranspiration and runoff, respectively.Plant must have sufficient woody biomass to support the mass of leaves, and the structure needs to be maintained by the proper allocation of carbon among stem (a s ), root (a R ), and leaf (a L ); the allocation of assimilated product is calculated using Equations ( 10)-( 12): where ε s , ε R , ε L and ω are plant function type (PET) dependent parameters following the rule of L and W represent light and water availability, respectively.

Developing GLOPEM Model and CEVSA Model : GLOPEM-CEVSA Model
GLOPEM-CEVSA represents a coupling of two independently developed models (GLOPEM and CEVSA), which describes the processes and mechanisms of an ecosystem based on remotely sensed vegetation indices (Figure 2, Table 2).In this model, the autophagic respiration simulation is based on plant physiological and ecological processes.The GLOPEM-CEVSA model reflects the vegetation growth process.By using the spatially continuous remote sensing data, the spatial heterogeneity of the ecosystem is considered in the model, which may improve the accuracy of NPP simulation [16].The most important feature of GLOPEM-CEVSA model is fPAR which could be retrieved by using remote sensing data to realize regional carbon flux simulation.In this paper, we improved the GLOPEM-CEVSA model, which uses MODIS atmosphere and land parameters to estimate model parameter fPAR and regional NPP.The model is based on a simplification of general radiation transfer equations [26].The NDVI and the NDVI-derived Simple Ratio (SR) were used to estimate fPAR, which integrated fPAR NDVI model and fPAR SR would provide an improved estimate [27].The vegetation classification accuracy is taken into account in the simulation of fPAR, which highlights the difference in photosynthetically active radiation absorption of different vegetation types.
Meanwhile, the differences among water vapor pressure, PAR, air temperature, the concentration of CO 2 , and soil moisture are taken into account in the light efficiency model [28,29].This approach can reflect the physiological and ecological processes of vegetation.
Among them, the vegetation PAR absorption obtained by using the climatology radiation calculation method was used for the estimation.fPAR refers to vegetation PAR absorption, which is used to describe the canopy structure and the function of the energy material exchange rate in the process of biophysical variables.Wang [18] uses the MODIS fPAR product to drive the model, but, in this paper, the precision of fPAR inversion with a space of 1 km and time scale of one month, was improved by Equations ( 16)- (19).
It has been found that fPAR was linearly related with NDVI and SR, previous study has shown that the NDVI-based fPAR model (fPAR NDVI ) and SR-based fPAR model (fPAR SR ) can overestimate or underestimate fPAR respectively, and then Los (1998) combines fPAR NDVI and fPAR SR , which would provide exact assessment result [30].
SR represents the ratio vegetation index.Studies show a good linear relationship among fPAR, NDVI, and SR.fPAR NDVI and fPAR SR represent the estimated fPAR, according to the SR and NDVI, respectively.NDVI i,max and NDVI i,min correspond to the maximum and minimum values of the vegetation types of NDVI, respectively, and correspond to a certain type of vegetation NDVI of 5% and 95% of the underside of percentile, respectively; the same situation applies to SR i,min and SR i,max .Among them, i represents different vegetation types, such us evergreen coniferous forest, evergreen broadleaved forest, deciduous coniferous forest, deciduous broad-leaved forest, shrub, grassland, farmland, urban-land, and bareland.fPAR min and fPAR max are unrelated to vegetation types.fPAR min value is 0.001, and fPAR max value is 0.95.Combined with the data of Chongqing, the maximum and the minimum of NDVI and SR for different vegetation types are calculated.The data involved in these calculations have been normalized and harmonized to the same spatial scale and coordinate system.LUE is the product of potential maximum light use efficiency and is also the key link in the GLOPEM-CEVSA model.We estimated LUE using critical environmental factors known to impact stomatal conductance, according to Equation (20): where PLUE is the potential maximum light utilization, and the model is based on the actual situation in Chongqing.The PLUE is estimated by different types of vegetation (Table 3), and the value is based on the simulation results of BIOME-BGG model by Running [31] and the actual situation in the southwest region.S(Ta), S(VPD), S(SWC), S(PAR) and S(CO 2 ), respectively, represent the effects of air temperature (Ta), vapor pressure (VPD), soil moisture (SWC), solar radiation (PAR) and carbon dioxide (CO 2 ) on stomatal conductance of vegetation.The temperature effect (T) can be divided into two factors.
where Te represent the average temperature ( • C), and T opt represent the most suitable average temperature in the month ( • C).
where VPD is vapor pressure deficit (g•kg −1 ), Qw(T) is the saturated specific humidity at the air temperature, and q is the air relative humidity.
where EET(x, t) is the actual evapotranspiration, and PET(x, t) is the potential evapotranspiration; PET(x, t) was provided by Holdridge [32].The actual EET is estimated by using the regional real evapotranspiration model [33].The solar radiation factor S(PAR) was calculated according to photosynthetic effective radiation by vegetation.
where CO 2 is the concentration of carbon dioxide.Usually, the autotrophic respiration (Ra) is divided into sustaining respiration (R m ) and growth respiration (R g ).
R g = r g GPP (30) where i represents different plant organs, i = 1, 2, 3, 4, representing roots, fine roots, stems, and leaves, respectively.M i is the biomass of each organ of vegetation, expressed as the concentration of C, r m,i for the vegetation of various organs of the maintenance of respiratory rate.Q 10 is the susceptibility factor of maintenance respiration to temperature, which means the change of maintenance respiration rate at 10 • C. In the BEPS model, the Q 10 of leaf, branch and root respiration were 2.1, 1.5, and 1.9, respectively.T is the temperature.Tb is the base temperature, using the BEPS model given T b = 15 • C. r g is the total growth of plant respiration in the proportion of total growth coefficient (r g = 0.35).
The autoregressive respiration of the original GLOPEM model that does not differentiate organs or contain biomass allocation simulation is modeled as an empirical relationship between vegetation biomass and air temperature [18,34].This process based autotrophic respiration simulation, better reflects the vegetation growth process, which may improve the accuracy of NPP simulation.

NPP Spatial Distribution and Change
This work adopted a new method to estimate regional NPP by using MODIS atmospheric and land parameters.This method improved the APAR and NPP calculation accuracy by using remote sensing parameters instead of coefficient transfer.The NPP from 2001 to 2011 was simulated by the remote sensing-process coupling model (GLOPEM-CEVSA).The average annual NPP distribution for 2001-2011 was plotted (Figure 3).According to the results of model inversion, we can find that the total annual NPP in the study area ranged from 0 to 1040.3 gC/m 2 a, the annual total NPP was nearly 0.5 PgC a −1 , and the mean NPP was 534.The overall regional changes indicated a high NPP in the east and the north and a low NPP in the west and the south.The spatial distribution of simulated NPP was corresponded to land cover type.The Daba Mountains and Wuling Mountains in the northeastern area in the central part of Chongqing are covered by extensive temperate forests, thereby resulting in a higher NPP that averaged more than 800 gC/m 2 a mostly in the Three Gorges reservoir area.This area has undergone extensive human disturbance, whereas the mountains are less disturbed as well as having a naturally more dense vegetation coverage.The western region and most of the main city were mainly low-value areas.The NPP between 100 and 200 gC/m 2 a is mainly along the coast of Yangtze River and urban built-up zone.Meanwhile, the NPP within the range of 300-700 gC/m 2 a was distributed in the central and western region in Chongqing.This is related to the distribution of karst topography in Chongqing.The overall regional changes indicated a high NPP in the east and the north and a low NPP in the west and the south.The spatial distribution of simulated NPP was corresponded to land cover type.The Daba Mountains and Wuling Mountains in the northeastern area in the central part of Chongqing are covered by extensive temperate forests, thereby resulting in a higher NPP that averaged more than 800 gC/m 2 a mostly in the Three Gorges reservoir area.This area has undergone extensive human disturbance, whereas the mountains are less disturbed as well as having a naturally more dense vegetation coverage.The western region and most of the main city were mainly low-value areas.The NPP between 100 and 200 gC/m 2 a is mainly along the coast of Yangtze River and urban built-up zone.Meanwhile, the NPP within the range of 300-700 gC/m 2 a was distributed in the central and western region in Chongqing.This is related to the distribution of karst topography in Chongqing.

The Distribution in Different Vegetation Types
From 2001 to 2011, NPP increased in shrub, evergreen broadleaved forest, and farmland (Table 4), whereas it decreased in deciduous coniferous forest, evergreen coniferous forest, coniferous forest vegetation, and bare land.The leaves of broadleaved forest areas are large; therefore, PAR absorption is strong, which resulted in the highest average NPP per unit area.The annual mean NPP of shrub was higher than that of coniferous forest, which may be explained by the low leaf area index of coniferous forest, the small leaf area, and the low utilization degree of light and heat, whereas the shrub area was luxuriant and therefore had strong absorption.The average NPP of grassland and cropland was 300 gC/m 2 a.However, the contribution of NPP of cropland to the study was higher because of the relatively large area of cropland.

The Distribution in Different Vegetation Types
From 2001 to 2011, NPP increased in shrub, evergreen broadleaved forest, and farmland (Table 4), whereas it decreased in deciduous coniferous forest, evergreen coniferous forest, coniferous forest vegetation, and bare land.The leaves of broadleaved forest areas are large; therefore, PAR absorption is strong, which resulted in the highest average NPP per unit area.The annual mean NPP of shrub was higher than that of coniferous forest, which may be explained by the low leaf area index of coniferous forest, the small leaf area, and the low utilization degree of light and heat, whereas the shrub area was luxuriant and therefore had strong absorption.The average NPP of grassland and cropland was 300 gC/m 2 a.However, the contribution of NPP of cropland to the study was higher because of the relatively large area of cropland.

Monthly Variation of NPP Value
Monthly distributions of NPP reflect interactions between physical (e.g., FPAR, precipitation) and biological (e.g., species composition, microbial activity, and ecological interactions) processes.The total NPP in the study area accumulated mainly from May to September (accounting for 84.2% of the annual NPP) and peaked in July and August (accounting for 44.3% of t annual NPP) (Figure 5).The following are three possible explanations for the seasonal distributions: (1) rainfall was not a limiting factor for vegetation growth in summer; (2) the temperature was closer to optimum in summer than in other seasons; and (3) FPAR peaked in summer, which led to the greatest fraction of APAR.
Calculated by the ratio between the months NPP difference for 11 years, the monthly rate change of NPP was obtained.From the change rate of monthly NPP (Figure 6), large variations occurred from March to May, and this period is the fastest growing period of NPP in Chongqing, with significant changes, whereas there is little change in the summer, and stability after September.Among the 11 years, the change rate of NPP in March reached 62.28%, while the lowest is 8.28% in August.This phenomenon indicates that NPP is stable in summer and the monthly variation of NPP in spring is larger, and the result may be caused by meteorological conditions in different seasons.total NPP in the study area accumulated mainly from May to September (accounting for 84.2% of the annual NPP) and peaked in July and August (accounting for 44.3% of t annual NPP) (Figure 5).The following are three possible explanations for the seasonal distributions: (1) rainfall was not a limiting factor for vegetation growth in summer; (2) the temperature was closer to optimum in summer than in other seasons; and (3) FPAR peaked in summer, which led to the greatest fraction of APAR.Calculated by the ratio between the months NPP difference for 11 years, the monthly rate change of NPP was obtained.From the change rate of monthly NPP (Figure 6), large variations occurred from March to May, and this period is the fastest growing period of NPP in Chongqing, with significant changes, whereas there is little change in the summer, and stability after September.Among the 11 years, the change rate of NPP in March reached 62.28%, while the lowest is 8.28% in August.This phenomenon indicates that NPP is stable in summer and the monthly variation of NPP in spring is larger, and the result may be caused by meteorological conditions in different seasons.

Spatial and Temporal Comparison among NPP Values
To validate the simulated NPP, two verification methods are used in the study: field observation, which is the most convincing method; and cross-validation between the relevant models, especially with the MODIS NPP product MOD17 correlation between the comparison.As we all know, measured NPP data for large areas are not available; however, vegetation productivity and biomass are strongly correlated with NPP [35].To validate the simulated NPP, the relationship  Calculated by the ratio between the months NPP difference for 11 years, the monthly rate change of NPP was obtained.From the change rate of monthly NPP (Figure 6), large variations occurred from March to May, and this period is the fastest growing period of NPP in Chongqing, with significant changes, whereas there is little change in the summer, and stability after September.Among the 11 years, the change rate of NPP in March reached 62.28%, while the lowest is 8.28% in August.This phenomenon indicates that NPP is stable in summer and the monthly variation of NPP in spring is larger, and the result may be caused by meteorological conditions in different seasons.

Spatial and Temporal Comparison among NPP Values
To validate the simulated NPP, two verification methods are used in the study: field observation, which is the most convincing method; and cross-validation between the relevant models, especially with the MODIS NPP product MOD17 correlation between the comparison.As we all know, measured NPP data for large areas are not available; however, vegetation productivity and biomass are strongly correlated with NPP [35].To validate the simulated NPP, the relationship

Spatial and Temporal Comparison among NPP Values
To validate the simulated NPP, two verification methods are used in the study: field observation, which is the most convincing method; and cross-validation between the relevant models, especially with the MODIS NPP product MOD17 correlation between the comparison.As we all know, measured NPP data for large areas are not available; however, vegetation productivity and biomass are strongly correlated with NPP [35].To validate the simulated NPP, the relationship between aboveground biomass and NPP can be applied.Biomass was measured in the field at 52 plots in September 2013 (Figure 7A).Field surveys were conducted on arbors (Figure 7B), shrubs (Figure 7C), and herbs (Figure 7D), with 23 arbor samples, 21 shrub samples, and eight herb samples.The plots were divided into four 5 m × 5 m plots, and the two plots of the diagonal were selected to investigate the shrub layer.The herbaceous community randomly set a 1 m × 1 m small plot in the grassland for the herbal investigation.The latitude and longitude of the survey plots and environmental factors such as slope, slope direction, soil type, and vegetation restoration time were recorded.The position of each field sample is measured by GPS, and the pixels of vegetation in the sample are pure when it is selected, so that the sample of the verification point and the pixel on the remote sensing image correspond.The ratio of belowground NPP (BNPP) to aboveground NPP (ANPP) was calculated using the following equations.
where BGB is totally accumulated by belowground biomass (BGB), nBGB refers to the BGB for a current year, and turnover is the vegetation root turnover [36].BGB and ANPP are provided by field measurements.We set liveBGB/BGB to 0.79 based on a field study by [37].BNPP/ANPP ratios were 9.75, 7.25, and 3.29 in arbor, shrub, and herb communities, respectively.Dry weight of biomass was converted to weight of carbon using carbon content.Carbon content was considered 45% for aboveground biomass [38,39].
where BGB is totally accumulated by belowground biomass (BGB), nBGB refers to the BGB for a current year, and turnover is the vegetation root turnover [36].BGB and ANPP are provided by field measurements.We set liveBGB/BGB to 0.79 based on a field study by [37].BNPP/ANPP ratios were 9.75, 7.25, and 3.29 in arbor, shrub, and herb communities, respectively.Dry weight of biomass was converted to weight of carbon using carbon content.Carbon content was considered 45% for aboveground biomass [38,39].
In a comparison with field data, the correlation between measured NPP and estimated NPP was 0.51 (Figure 8), thereby showing a significant linear relationship between measured NPP and simulated NPP (p < 0.01).
In a comparison with field data, the correlation between measured NPP and estimated NPP was 0.51 (Figure 8), thereby showing a significant linear relationship between measured NPP and simulated NPP (p < 0.01).This situation occurs because most of Chongqing's districts are located in the Three Gorges reservoir area, which has high forest coverage and is well-protected.The results of this study lie between those of the improved CASA model (R 2 = 0.756, Figure 9A) and the CEVSA model, which is close to the estimation result of Dong Dan, in which the vegetation characteristics of karst area are considered.The correlation between estimated NPP and MOD17 NPP was 0.56 (Figure 9B), indicating that the model implies a high accuracy.The spatial variation trend is consistent, but the NPP of evergreen coniferous forest was slightly higher.Miami and Thorthwait models are more accurate for simulating grass NPP, which is similar to the result of study.In this paper, the estimation of shrubs is not very good because of the extensive distribution of karst areas in Chongqing and the carbonate exposed area, which covered 32,038.14km 2 and accounted for 38.9% of the total area.The corrosion and erosion of the rock are serious, which affects the development of low trees.In sum, the results are in a reasonable range, and the spatial variation trend is consistent because of the difference in time scale, spatial region scale, and vegetation type.This situation occurs because most of Chongqing's districts are located in the Three Gorges reservoir area, which has high forest coverage and is well-protected.The results of this study lie between those of the improved CASA model (R 2 = 0.756, Figure 9A) and the CEVSA model, which is close to the estimation result of Dong Dan, in which the vegetation characteristics of karst area are considered.The correlation between estimated NPP and MOD17 NPP was 0.56 (Figure 9B), indicating that the model implies a high accuracy.The spatial variation trend is consistent, but the NPP of evergreen coniferous forest was slightly higher.Miami and Thorthwait models are more accurate for simulating grass NPP, which is similar to the result of study.In this paper, the estimation of shrubs is not very good because of the extensive distribution of karst areas in Chongqing and the carbonate exposed area, which covered 32,038.14km 2 and accounted for 38.9% of the total area.The corrosion and erosion of the rock are serious, which affects the development of low trees.In sum, the results are in a reasonable range, and the spatial variation trend is consistent because of the difference in time scale, spatial region scale, and vegetation type.
accurate for simulating grass NPP, which is similar to the result of study.In this paper, the estimation of shrubs is not very good because of the extensive distribution of karst areas in Chongqing and the carbonate exposed area, which covered 32,038.14km 2 and accounted for 38.9% of the total area.The corrosion and erosion of the rock are serious, which affects the development of low trees.In sum, the results are in a reasonable range, and the spatial variation trend is consistent because of the difference in time scale, spatial region scale, and vegetation type.

Relationships between NPP and Meteorological Factors
The relationships among NPP, temperature and precipitation have been previously reported.Results showed that NPP was correlated with temperature (R 2 = 0.924, Figure 10A) and precipitation (R 2 = 0.571, Figure 10B) at every month interval.The monthly change in NPP was similar to that in precipitation, especially from January to June and from August to December.NPP continued to increase in July, although precipitation was lower than that in July (Figure 10B), whereas NPP and temperature had the same trend (Figure 10A).The highest temperature occurred in July, during which NPP also peaked.The probable reason for this situation was that water-absorbing processes caused a lag in vegetation growth in the study area (Figure 11).In this paper, we analyzed the relationship between NPP and meteorological factors without considering vegetation cover.We can find that NPP was consistent with the change of temperature and was not affected by precipitation fluctuation, which indicated that the value of NPP was mainly affected by temperature.NPP sharply decreased from July to August.The spatial distribution of NPP could also be explained by the relationship to temperature/precipitation via the spatial distribution of heat and water in the Chongqing area (Figure 11), and the special mountainous terrain conditions have a significant effect on the distribution of NPP.

Relationships between NPP and Meteorological Factors
The relationships among NPP, temperature and precipitation have been previously reported.Results showed that NPP was correlated with temperature (R 2 = 0.924, Figure 10A) and precipitation (R 2 = 0.571, Figure 10B) at every month interval.The monthly change in NPP was similar to that in precipitation, especially from January to June and from August to December.NPP continued to increase in July, although precipitation was lower than that in July (Figure 10B), whereas NPP and temperature had the same trend (Figure 10A).The highest temperature occurred in July, during which NPP also peaked.The probable reason for this situation was that water-absorbing processes caused a lag in vegetation growth in the study area (Figure 11).In this paper, we analyzed the relationship between NPP and meteorological factors without considering vegetation cover.We can find that NPP was consistent with the change of temperature and was not affected by precipitation fluctuation, which indicated that the value of NPP was mainly affected by temperature.NPP sharply decreased from July to August.The spatial distribution of NPP could also be explained by the relationship to temperature/precipitation via the spatial distribution of heat and water in the Chongqing area (Figure 11), and the special mountainous terrain conditions have a significant effect on the distribution of NPP.A B   FPAR was another factor that influenced spatial and temporal variation in NPP.Different land cover types varied in FPAR patterns, thereby generating a variation in PAR absorption and consequently changing NPP.Moreover, seasonal changes in FPAR in the IMAR caused temporal NPP changes.In this study, FPAR inversion in particular was based on the satellite remote sensing data.Thus, the relationship between FPAR and atmospheric parameters and vegetation mechanism can be seen more easily when using MODIS than when only meteorological station statistical data are used.

Sensitivity Analysis
The sensitivity coefficient was the rate of change of the ratio between NPP and meteorological factors, such as precipitation, temperature, and sunshine hours.
S V i , which is non-dimensional, is the sensitivity coefficient of V i .NPP and ∆NPP are the NPP and the rate of change, respectively.V i and ∆V i are meteorological factors and the rate of change, respectively.A greater sensitivity coefficient corresponds to greater influence of the meteorological factors on NPP.The annual variation of sensitivity coefficient (Figure 12) also has unique distribution characteristics.Judging from the characteristics of multi-year changes, the sensitivity coefficient of the temperature and precipitation is dominated by stable fluctuation.The annual changes are not obvious, and the sensitivity coefficient of the temperature fluctuates around 0. NPP changes.In this study, FPAR inversion in particular was based on the satellite remote sensing data.Thus, the relationship between FPAR and atmospheric parameters and vegetation mechanism can be seen more easily when using MODIS than when only meteorological station statistical data are used.

Sensitivity Analysis
The sensitivity coefficient was the rate of change of the ratio between NPP and meteorological factors, such as precipitation, temperature, and sunshine hours.
S V i , which is non-dimensional, is the sensitivity coefficient of V i .NPP and ∆NPP are the NPP and the rate of change, respectively.V i and ∆V i are meteorological factors and the rate of change, respectively.A greater sensitivity coefficient corresponds to greater influence of the meteorological factors on NPP.The annual variation of sensitivity coefficient (Figure 12) also has unique distribution characteristics.Judging from the characteristics of multi-year changes, the sensitivity coefficient of the temperature and precipitation is dominated by stable fluctuation.The annual changes are not obvious, and the sensitivity coefficient of the temperature fluctuates around 0.63.The sensitivity coefficient of precipitation decreased from 2007 (the sensitivity coefficient of 2007 is 0.48).The sensitivity coefficient of sunshine hours showed an upward trend.The fluctuation is small from 2003 to 2005, and the sensitivity coefficient is approximately 0.4.This finding indicates that the sensitivity of the sunshine increased over the 11 years.

Uncertainties of the Model and the Interference of Mixed Pixels
In the study, based on MODIS NDVI data, land parameters and meteorological data derived from 34 meteorological stations and from Chongqing for years 2001-2011, we improved the GLOPEM-CEVSA model and investigated the seasonal variation and spatial distribution of NPP as well as analyzed the relationship between NPP and the climate factors.However, our estimate still retain some uncertainties for the following reasons.
First, there exists the scale effect between remotely sensed data and meteorological data.The spatial resolution of MODIS NDVI data is 1 KM, affecting the inversion accuracy of FPAR.In the process of interpolating meteorological data, we have taken the average value to reduce the estimation error caused by the spatially or temporally scale effect.Second, the validated data used in the study came from the MOD17 NPP, which is a global land vegetation net productivity product calculated by model and has been validated widely, across most of the world [40].In the process of the calculations, some of the maintenance respiration could consume the output of the total NPP, which is termed PSNnet (net photosynthesis).Thus, it is necessary to distinguish the true NPP and forest inventory data.NPP, as a result of the verification is modeled by the measured biomass data, and the biomass data collected in the field inevitably has some error, resulting in some errors for the verification results.
Understanding and identifying the sources of uncertainties and then devoting efforts to improving them are critical in improving the accuracy of model estimation; therefore, more research is needed in the future for reducing the uncertainties from different sources in the process of NPP simulation, especially for the localization of model parameters, spatial scale matching between each module data of the model.
The mixed pixel is often one of the most important error sources for the simulation of NPP.Because the land cover mixtures are more complex and spatially heterogeneous, most pixels of MODIS data consist of several land cover due to the low spatial resolution, which will complicate the quantitative evaluation of mixed pixel effect on NPP, which affects the estimation accuracy of net primary productivity to a certain extent.Therefore, NPP assessment with complex land cover mixtures still needs further investigation.

Conclusions
This study took Chongqing as a case study, and used meteorological data, and MODIS atmospheric and land parameters to develop a GLOPEM-CEVSA model to estimate NPP (net primary productivity), and then analyzed the spatial and temporal distribution of NPP for different vegetation types.The NPP from April to October accounted for 80.1% of the annual NPP, while that from June to August accounted for 43.2%.The annual NPP in the study area was approximately 534 gC/m 2 a from 2001 to 2011, and showed a gradually decreasing trend from northeast to southwest.The average NPP of forest, scrub, cropland, grassland, urban land, and bare land was 412.29, 554.83, 347.93, 362.40, 203.32 and 110.68 gC/m 2 a, respectively, and the NPP showed differences among different land cover types.
The GLOPEM-CEVSA model not only improved the accuracy of PAR and NPP estimates by using remote sensing parameters instead of coefficient transfer, but also integrated MODIS land parameters to drive the model.The growth rate of vegetation was mainly determined by precipitation and temperature factors.NPP in the karst region was significantly correlated with precipitation (R 2 = 0.571) and temperature (R 2 = 0.924) at monthly temporal scales and showed stronger sensitivity to interannual variation in temperature.Although the study showed promising results for remote sensing based vegetation NPP estimation, there are limitations to the accuracy of fPAR inversion using MODIS data; therefore, further work is needed to improve the estimation accuracy.In addition, spatio-temporal patterns of NPP and their relationships with climate factors still need further research in the Karst area.Chen et al. [41] have already studied the relationship among precipitation, relative humidity, mean temperature and vegetation EVI in Karst rocky area.However, there are many other meteorological factors, such as how the growing season maximum and minimum temperatures affect NPP.Therefore, the vegetation response to climate change is complex, and deserves more detailed and deeper inquiry in future research.

Figure 1 .
Figure 1.Meteorological location and topography of Chongqing.

Figure 1 .
Figure 1.Meteorological location and topography of Chongqing.

Figure 3 .
Figure 3. Spatial distribution of mean NPP from the year 2001 to 2011.

Figure 3 .
Figure 3. Spatial distribution of mean NPP from the year 2001 to 2011.
The NPP value appeared to fluctuate from the year 2001 to 2011 in the whole area (Figure 4), the range of NPP is 465.60~646.11gC/m 2 a, with an average of 534.55 gC/m 2 a.In 2003 and 2008, the values show a significant peak, with annual NPP average of 646.11 gC/m 2 a and 555.98 gC/m 2 a. Especially the NPP value of 2003 is the highest among the 11 years.The values of NPP is relatively low for the years 2001, 2006 and 2009, the mean values of NPP are 465.60 gC/m 2 a, 504.45 gC/m 2 a and 517.47 gC/m 2 a, respectively.The total variation range of the NPP value is 3.83 × 10 13 ~5.32 × 10 13 gC/a, the average is 4.4 × 10 13 gC/a, and the fluctuation range is 0.59 × 10 13 gC/a, accounting for 13.4% of the total.

Figure 4 .
Figure 4.The NPP value from the year 2001 to 2011.

Figure 4 .
Figure 4.The NPP value from the year 2001 to 2011.

Figure 6 .
Figure 6.The change rate of monthly NPP from the year 2001 to 2011.

Figure 6 .
Figure 6.The change rate of monthly NPP from the year 2001 to 2011.

Figure 6 .
Figure 6.The change rate of monthly NPP from the year 2001 to 2011.

Figure 8 .
Figure 8.Comparison of the values between the measured NPP and estimated NPP.

Figure 9 .Figure 8 .
Figure 9.Comparison of NPP values among the estimated model, the CASA model (A) and MOD17

Figure 9 .Figure 9 .
Figure 9.Comparison of NPP values among the estimated model, the CASA model (A) and MOD17 (B).

Figure 10 .
Figure 10.Relationships between estimated NPP and precipitation and temperature (A,B).

Figure 10 .
Figure 10.Relationships between estimated NPP and precipitation and temperature (A,B).

Figure 10 .
Relationships between estimated NPP and precipitation and temperature (A,B).

Figure 11 .Figure 11 .
Figure 11.Monthly changes of estimated NPP and precipitation (A), temperature (B) (from 2001 to 2011).FPAR was another factor that influenced spatial and temporal variation in NPP.Different land cover types varied in FPAR patterns, thereby generating a variation in PAR absorption and consequently changing NPP.Moreover, seasonal changes in FPAR in the IMAR caused temporal 63.The sensitivity coefficient of precipitation decreased from 2007 (the sensitivity coefficient of 2007 is 0.48).The sensitivity coefficient of sunshine hours showed an upward trend.The fluctuation is small from 2003 to 2005, and the sensitivity coefficient is approximately 0.4.This finding indicates that the sensitivity of the sunshine increased over the 11 years.Remote Sens. 2017, 9, 477 17 of 22

5. 3 .
Uncertainties of the Model and the Interference of Mixed Pixels In the study, based on MODIS NDVI data, land parameters and meteorological data derived from 34 meteorological stations and from Chongqing for years 2001-2011, we improved the GLOPEM-CEVSA model and investigated the seasonal variation and spatial distribution of NPP as

Table 2 .
The key conceptions and parameters in the framework of model.
NPP Net primary productivity (gC/m 2 a) PAR Photosynthetically active radiation (MJ/m 2 a) fPAR Photosynthetically active radiation LUE Actual light use efficiency PLUE Potential maximum light use efficiency Ra Autotrophic respiration S(PAR) Radiation stress S(T) Temperature stress S(SWC) Soil water stress S(VPD)Vapor stress

Table 2 .
The key conceptions and parameters in the framework of model.

Table 3 .
The PLUE of different LUCC.

Table 4 .
The statistics of average NPP in different vegetation types.