Optimizing the photosynthetic parameter Vcmax by assimilating MODIS-fPAR and MODIS-NDVI with a process-based ecosystem model

Agricultural and Forest Meteorology 198–199 (2014) 320–334

Contents lists available at ScienceDirect

Agricultural and Forest Meteorology journal homepage: www.elsevier.com/locate/agrformet

Optimizing the photosynthetic parameter Vcmax by assimilating MODIS-fPAR and MODIS-NDVI with a process-based ecosystem model Shi Hu 1 , Xingguo Mo ∗ , Zhonghui Lin 2 Key Laboratory of Water Cycle and Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, PR China

a r t i c l e

i n f o

Article history: Received 9 September 2013 Received in revised form 30 August 2014 Accepted 7 September 2014 Keywords: Vcmax Interval estimation North China Plain MODIS VIP model

a b s t r a c t Combination of satellite remote sensing data and an ecosystem model provides an opportunity to monitor net ecosystem production, water cycle and energy balance at the regional scale. Photosynthesis is a critical ecological process that is coupled to the carbon and water cycle and energy balance. Therefore, an accurate description of its spatiotemporal pattern is essential when simulating an ecosystem at the regional scale. To determine the spatial distribution of the maximum Rubisco catalytic capacity (Vcmax ), we have developed a scheme that optimizes the photosynthetic parameter from a remotely sensed fPAR (the fraction of photosynthetically active radiation absorbed by the plant canopy) and NDVI (normalized difference vegetation index) using the VIP ecosystem model. It integrates the interval estimation method and a one-dimensional searching algorithm, in which the samples include randomly selected pixels, the photosynthetic capabilities are optimized with the golden section search algorithm in the randomly sampled pixels to derive the prior probability of Vcmax , and then the search interval of Vcmax is narrowed to a confidence interval. We verified this scheme on the North China Plain (NCP) to determine the Vcmax pattern in winter wheat at a 1-km resolution. The simulation results by the VIP model with the derived Vcmax pattern were indirectly validated using census data for grain yield, field evapotranspiration (ET) measurements, the MODIS leaf area index (MODIS-LAI) and daily MODIS land surface temperatures (MODIS-LST). The validation results demonstrated a satisfactory agreement between the simulated and measured data with R2 of 0.63, 0.82, 0.29 and 0.92 for yield, ET, LAI and LST, respectively. It is suggested that the proposed Vcmax -retrieving method is practical for regional crop growth predictions. © 2014 Elsevier B.V. All rights reserved.

1. Introduction The combination of satellite remote sensing and ecosystem models provides an opportunity for monitoring net ecosystem production, the water cycle, and energy balance on a regional level (Dorigo et al., 2007). Photosynthesis is a critical ecological process that is coupled to the ecosystem’s carbon cycle, water cycle and energy balance. Therefore, an accurate description of its spatiotemporal pattern is essential when simulating an ecosystem at the regional scale. Biophysical processes, which include CO2 transport through the leaf and stomata, and biochemical processes in the chloroplast thylakoid membranes, stroma, mitochondria and cytosol of the cell

∗ Corresponding author. Tel.: +86 10 64889307; fax: +86 10 64851844. E-mail addresses: [email protected] (S. Hu), [email protected] (X. Mo), [email protected] (Z. Lin). 1 Tel.: +86 10 64888101. 2 Tel.: +86 10 64889063. http://dx.doi.org/10.1016/j.agrformet.2014.09.002 0168-1923/© 2014 Elsevier B.V. All rights reserved.

determine the net rate of CO2 assimilation (A). These biophysical and biochemical processes in conjunction with environmental variables, such as light intensity and temperature, confer different effects on A; thus, it is difficult to predict how A will be affected by genetics, epigenetics and the environment. Therefore, certain canopy-state variables, such as the leaf area index (Doraiswamy et al., 2004; Moulin et al., 2003) and plant chlorophyll concentration (Haboudane et al., 2002; Zhao et al., 2004), which are closely related to the plant photosynthetic capability, are frequently used to describe A. The canopy state variables can be easily retrieved from remote sensing data using the canopy reflectance model (Shabanov et al., 2000; Schlerf and Atzberger, 2006) or the statistical relationship between the spectral signature and biophysical (or biochemical) variables (Atzberger et al., 2003; Huang et al., 2004). Thus, the retrieved state variables are directly used to replace variables in the model because the observed data are presumably more reliable than simulated data. This approach is widely employed in studies that range from sub-regional (Bastiaanssen and Ali, 2003) to regional (or continental) scales (Mo et al., 2005) due to its excellent performance and its simplicity. However, the observed errors in the

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

remote sensing state variables are propagated into the model when the remote sensing state variable is used directly (Dorigo et al., 2007). This shortcoming can be avoided by adjusting the model parameters to determine the optimal consistency between the simulated and observed remote sensing state variables (Launay and Guerif, 2005). Parameter adjusting is more flexible for assimilating remote sensing state variables and their associated error in the model because a physical description of the underlying process is an acceptable representation of the natural system, and the adjusted model parameter is expected to improve model prediction. Therefore, using remote sensing data to adjust the model parameters has been the more frequent approach instead of direct use of remote sensing state variables in a model. This practice was confirmed by the following groups. Dente et al. (2008) used remotely sensed LAI to improve wheat yield simulation and determine an optimal set of input parameters, including sowing data, soil wilting point and field capacity. Fang et al. (2008) utilized the MODIS-LAI to reinitialize planting data, planting density and row spacing as well as improve the model’s performance in simulating the corn yield. Hazarika et al. (2005) successfully simulated net primary productivity using the re-initialized leaf area index. However, no reports have described photosynthetic parameters for re-initialization due to the complexity of the photosynthetic process. The most sensitive photosynthetic parameters in the photosynthesis process are Vcmax (the maximum rate of Rubisco carboxylase activity) and Jmax (the maximum rate of photosynthetic electron transport). Vcmax describes a photosynthesis limitation, the photosynthetic Rubisco enzyme system, which dominates under light saturation as regulated by the photosynthetically active radiation (PAR) effectively absorbed by photosystem II (von Caemmerer and Farquhar, 1981). Vcmax and Jmax are tightly coupled, and the fixed ratio of Jmax /Vcmax at 25 ◦ C is typically assumed in largescale modeling schemes. Therefore, specifying Vcmax is critical and substantially contributes to terrestrial GPP (gross primary productivity) sensitivity to the atmospheric CO2 concentration and the overall uncertainties of model-predicted carbon fluxes. In wheat, a 40% decrease in Vcmax led to a 34% decrease in GPP (Lei et al., 2011), and a 10% increase in Vcmax led to a 4% increase in yield in NCP (Mo et al., 2005). Although both parameters are tightly correlated with a plant’s photosynthetic capacity, carboxylation capacity (Vcmax ) has been measured and studied more extensively (Wullschleger, 1993; Leuning, 1997; Kattge and Knorr, 2007). At the leaf scale, the key photosynthetic parameters, Vcmax and Jmax , can be calculated by fitting the CO2 response curve to the measured photosynthetic rate and intercellular CO2 concentration (Muller et al., 2005; Alonso et al., 2009). At the regional scale, the photosynthetic rate and intercellular CO2 concentration are unavailable for each pixel; thus, it is impossible to obtain photosynthetic parameter patterns using a CO2 response curve. Vcmax is affected by factors that determine the level of activated enzyme (Rubisco), which is determined by the foliar chlorophyll concentrations. These factors include temperature (Sage, 2002), water stress (Flexas et al., 2006), leaf nitrogen content (Maroco et al., 2002; Gonzalez-Real and Baille, 2000), radiation intensity (Choi and Roh, 2003), salt stress (Centritto et al., 2003) and a combination of the aforementioned factors. Under environmental stress, chlorophyll levels decrease and induce a decrease in Vcmax (Houborg et al., 2013), because leaf chlorophyll levels are the primary factor that control spectral reflectance variation in the visible light regions and dominantly control the solar radiation levels that a leaf absorbs. The decrease in chlorophyll levels affects the abilities of green vegetation to reflect incident solar radiation; this response can be detected using the remote sensing data, such as: MODIS-fPAR and MODIS-NDVI (Dawson et al., 2003; Huang et al., 2010). Therefore, stressed vegetation or vegetation with a low photosynthetic capability has a reduced fPAR value (or

321

NDVI value), which suggests that some remote sensing data are reliable indicators of photosynthetic capability in specific vegetation. Therefore, developing an empirical or theoretical transfer function between Vcmax and remote sensing data (such as: MODIS-fPAR or MODIS-NDVI) could be an approach for obtaining spatial patterns for photosynthetic parameters. The quantitative relationship between Vcmax and remote sensing data recently has been the subject of interest. Many of the conditions that affect plant development and grain yield either favorably or adversely (i.e. fertilization treatment, rust infection, drought or precipitation events) result in a corresponding increase or reduction of the crop’s photosynthetically active biomass at the canopy scale, and this response can often be captured though spectral measures such as the normalized difference vegetation index (NDVI) (Tucker, 1979). However, the NDVI often fails to reflect plant photosynthetic capacity accurately, especially when the time of NDVI collection does not coincide with the flourishing stage of the crop’s photosynthetically active biomass. Most authorities recommend collection of spectral data when the crop canopy is full (Basnyat et al., 2004; Becker-Reshef et al., 2010), because image acquisition before full crop canopy coverage results in spectral information that includes the contribution from the soil surface rather than from just the crop itself. Usually, the cumulative NDVI in the flourishing stage is a good indicator of plant photosynthetic capability (Basnyat et al., 2004), since it provides the aggregated assessment of plant photosynthetic capacity over time. Hu and Mo (2012) derived a regression between the accumulated NDVI during the growing season and the county-level statistic grain yield, and then Vcmax estimated using the process-based ecosystem model VIP with the assumption that the model has good performance in yield simulation at regional scale. This practice is similar to a statistical approach, it is restricted by the census data quality as well as heterogeneity of the soil and vegetation at the county level. In the data sparse region which lack of grain census data, the relationship between remote sensing data and photosynthetic capacity is hard to establish, and limited the use of the statistical approach. This study aims to develop a new scheme to retrieve the photosynthetic parameter (Vcmax ) for winter wheat over the NCP by assimilating MODIS-fPAR and MODIS-NDVI data with the VIP ecosystem model. This approach is trying to resolve the photosynthetic capacity pattern retrieval problem in data sparse region. The model predictions with the retrieved Vcmax pattern were evaluated using the census data for grain yield, measured ET, MODIS-LST and MODIS-LAI. Finally, the simulation error of the model with the retrieved Vcmax pattern and the difference between this scheme and other Vcmax retrieval approaches are discussed.

2. Methods 2.1. Model description The VIP model used in this study is a dynamic ecosystem model that simulates radiation, water, heat and CO2 transfer processes during the crop-growing season (Mo et al., 2012; Mo and Liu, 2001). The water cycle incorporates the total above-canopy evapotranspiration and soil moisture transfer. Energy fluxes are described with a two-source scheme that individually discerns the canopy and the soil surface. The temperatures of the canopy (Tc ) and soil surface (Tg ) which were used to calculate land surface temperature (Appendix A, Eq. (A3)) are governed by the twosources soil-canopy energy balance sub-model (Shuttleworth and Wallace, 1985; Mo et al., 2012). The carbon cycle includes assimilation via photosynthesis, crop growth and soil organic matter decomposition schemes. The leaf photosynthesis rate is primarily restricted by the carboxylation rate, which is limited by the

322

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

enzymatic activity of Rubisco, the photoelectron transfer rate, and the phosphate release rate during the consumption of photosynthetic products. The actual photosynthesis rate depends on the limit values of the three aforementioned aspects (Farquhar et al., 1980; Sellers et al., 1992). The net assimilation rate of CO2 (An , ␮mol C m−2 s−1 ) is therefore expressed as follows: An = min(AR , AE , AS ) − Rd

(1)

where AR is the Rubisco limiting rate, AE is the light limiting rate, AS is the capacity for the export or utilization of the photosynthetic products, and Rd is the leaf respiration rate (␮mol C m−2 s−1 ). The Rubisco limiting rate, AR , is expressed as follows (Farquhar et al., 1980; Sellers et al., 1992): AR = Vc



ci −  ∗ ci + Kc (1 + O2 /Ko )



(2)

Vc is Rubisco’s maximum catalytic capacity where (␮mol C m−2 s−1 ), ci and O2 are the intercellular concentrations of CO2 and O2 (Pa), respectively,  * is the CO2 compensation point (Pa, =0.5O2 /S, where S is the Rubisco specificity for CO2 relative to O2 ), and Kc and Ko are the Michaelis constants for CO2 and O2 , respectively (Pa). For more details regarding AE and AS , refer to Mo and Liu (2001). To account for the light extinction in a canopy, a multilayer scheme for both sunlit and shaded groups is developed to upscale the leaf photosynthesis to the entire canopy (Mo and Liu, 2001). In the VIP model, the canopy is divided into 20 layers for the photosynthesis calculations with the same leaf area index. In the sub-layer i, the Rubisco’s maximum catalytic capacity–Vc(i) is linearly related to the leaf nitrogen concentration that exponentially decreases in a canopy, therefore, Vc(i) should also decline exponentially with an extinction coefficient, kv , and is expressed as follows:



Vc(i) = Vcmax exp −kv LAI ×

i 20



(3)

where Vcmax is the value of Vc at the top of the canopy (i = 1), kv is taken as 0.6 and LAI is the leaf area index that is determined by the leaf dry mass in the dynamic crop growth process (Appendix B). By integrating the assimilation rate of all layers, the entire canopy assimilation rate is obtained. As stated above, the maximum catalytic capacity of Rubisco at the top layer of the canopy (Vcmax ) is a key parameter in the VIP model and derived from remote sensing data in this study. The value of Vcmax is affected by any factor that may affect the amount of activated enzyme (Rubisco) present, such as leaf age (Wilson et al., 2001), nitrogen starvation (Xu and Baldocchi, 2003) and drought (Misson et al., 2010), leading to spatial and temporal variation of Vcmax . During leaf senescence, Rubisco decomposes, and its concentration in the leaf falls, resulting in the decline and seasonal variation of Vcmax . As Vcmax is the mean growing season value in the VIP model, we would not pay more attention to temporal variation of Vcmax in this study. Nitrogen is another key factor in determining the photosynthetic capability and Vcmax value. The VIP model has not yet addressed the nitrogen cycle mechanically in this study, the influence of soil nitrogen is here assumed to be accounted for by Vcmax values. Thus, to some extent, Vcmax is considered to be an empirical parameter reflecting the background of crop photosynthetic capacity in this study.

on the simple and reliable fPAR merit function root mean square error (RMSE, Eq. (4)). The model produces the daily fPAR and photosynthetic parameter (Vcmax ) values when a minimization threshold is satisfied. However, the MODIS-fPAR is not available with a daily temporal resolution. The maximum simulated fPAR of eight days was retrieved to match the temporal resolution of the MODIS-fPAR . In Eq. (4), xi obs is the MODIS-fPAR , xi sim is the simulated fPAR , and n is the number of MODIS-fPAR observations.



RMSE =

n (x i=1 i obs

− xi

sim )

2

(4)

n

The golden section search algorithm (Press et al., 2007) used in this study was aimed at finding the minimum cost function by successively narrowing the Vcmax range wherein the minimum exists. Because the initial searching interval [40 ␮mol C m−2 s−1 , 120 ␮mol C m−2 s−1 ] was determined based on the special attributes of winter wheat, the assimilation strategy described above requires considerable computational time to obtain the optimal parameter value, especially for a large region. The Vcmax value is affected by factors that determine the level of activated enzyme (Rubisco); crops grown in similar environments typically have similar Vcmax values. Therefore, a Vcmax value from a similar environment defines a more narrow location interval than the value determined from its special attributes. The interval estimation, wherein sample data are used to calculate an interval of possible values for an unknown population parameter, provides a solution to identify a narrower search interval for the parameter. The interval containing a population parameter is established by applying knowledge on the fidelity of properties from a sample that represents the entire population. In other words, the interval that most likely includes the population parameter can be calculated using sample data from similar, corresponding environments. Based on the above rationale, a three-step technique was developed to assimilate MODIS-fPAR and the ecosystem model on a regional scale as follows: 1. Retrieval Vcmax value of sample. The sample contains 5% of the pixels (randomly selected) from the study area. Based on the merit function (RMSE) constructed with the simulated fPAR and MODIS-fPAR , the Vcmax value of the sample is retrieved using the golden section search algorithm. This initial step determines the relationship between the MODIS data and Vcmax value of the sample. 2. Interval estimation. Interval estimation is performed based on the relationship between the sample’s Vcmax value and its corresponding MODIS data. This step is performed to deduce the confidence interval of the population mean and standard deviation (SD) as well as narrow the population search interval in the searching algorithm. According to the central limit theorem and law of large numbers, regardless of whether a population is normally distributed, when the sample is large, the sample means obey the distribution N(, x2¯ ), and the variable (n − 1) 2 /s2 has a 2 distribution with (n − 1) degrees of freedom (DeGroot and Schervish, 2012). Therefore, at the significance level ˛ = 0.05, the confidence level of population mean and population standard deviation can be expressed as Eqs. (5) and (6), respectively. P( − 1.96x¯ < x¯ <  + 1.96x¯ ) = 0.95



2.2. Assimilation strategy In this article, “assimilation” may represent any combination of a model and remotely sensed observations. A general outline of our assimilation strategy is provided in Fig. 1. fPAR simulated by the crop model is compared with MODIS-fPAR products from the optimization process. The crop model reinitiates the input parameters based

P

(n − 1) 2 2˛/2 (n − 1)

< s2 <

(n − 1) 2 21−˛/2 (n − 1)

(5)

= 0.95

(6)

When the sample is large, and the population variance is unknown, the population mean’s 95% confidence interval (L1 mean , L2 mean ) can be generated using the sample mean x¯ and

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

323

Initial searching interval

VIP model

New searching interval [min, max]

Simulated fPAR

Modis-fPAR

Merit function f(Vcmax) Golden section searching algorithm

NO

Minimum reached?

YES Optimal Vcmax Fig. 1. Flowchart for the retrieval of the Vcmax value.

standard deviation x¯ with Eq. (7), the confidence interval of the population SD at a 95% confidence level (L1 SD , L2 SD ) can be generated using the sample SD with Eq. (8). L1

mean

L1

=

SD

= x¯ −



1.96x¯ , √ n

√ n−1 2˛/2 (n − 1)

L2 ,

mean

L2

SD

1.96x¯ √ n √ n−1

= x¯ + =



(7)



(8)

21−˛/2 (n − 1)

where L1 mean and L2 mean are the upper and lower confidence limits of population mean, L1 SD and L2 SD are the upper and lower confidence limits of the population SD, n is the number of samples. The 2 value at (n − 1) degrees of freedom can be generated using the 2 distribution tables (DeGroot and Schervish, 2012). 3. Acquisition of the Vcmax pattern. The confidence interval of population mean and SD can be used to generate a narrowed population search interval. Based on the new search interval, the Vcmax value of any pixel in the study area can be obtained using the golden section search algorithm. 2.3. Validation method The ecosystem’s crop-growth processes, water balance and energy balance are affected by the crop photosynthetic capability. The census data on yield, measured ET, MODIS-LST and MODISLAI were used to validate the ecosystem model simulation results

(VIP model) for the Vcmax pattern, as observation data on the photosynthetic parameters were unavailable at the regional scale, this practice indirectly verified the Vcmax pattern. Because certain wheat fields in the study area were too small for adequate representation by the 1 km spatial resolution RS data, the Vcmax value represents the average value of the mixed pixel in a 1 × 1 km grid, including crop and bare soil. To reduce the errors introduced by the mismatches between the RS product and ground data in the validation, we considered the wheat area percentage in each square in a 1 × 1km grid for the validation. The field data used in the validation (including the above-ground biomass, ET and LAI) were collected from the grid with a wheat field percentage greater than 95%.

3. Study area and data preparation 3.1. Study area The NCP is the second largest plain in China, with an area of 33 × l04 km2 . The NCP extends from 32◦ 00 N to 40◦ 24 N and from 112◦ 48 E to 122◦ 45 E. The plain has a sub-humid warm temperate continental monsoon climate with a mean annual temperature range of 10.8–14.1 ◦ C and an annual precipitation of 500–1000 mm, which is distributed irregularly among the seasons. Correspondingly, the prevailing cropping system is a wheat–maize rotation, which represents the traditional cropping system in the NCP and accounts for 90% of the cereal sowing area (Fig. 2). Due to the

324

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

Fig. 2. Location of the study area (depicted in red in the panel a), (a) the ratio of farmland, the location of 20 agricultural meteorological stations and three agricultural experimental stations and (b) soil texture in the NCP and the location of national meteorological stations in and around NCP.

insufficient rainfall in the wheat growing period (October–May), most farmlands (approximately 70%) are typically irrigated several times to ensure high productivity levels. 3.2. Data preparation Three types of data were used, namely, geographical, meteorological and agricultural data. The geographical data were used as the basic land surface information for the crop and soil. The meteorological data were used as atmospheric forces that drive the VIP model. The agricultural data were used to validate the model and crop parameters. 3.2.1. Land surface characterization data The land surface characterization includes the topography, soil physical texture map and land use map. The digital elevation model was obtained from the topographic map at a 1:250,000 scale, which can be created using a terrain raster with a resolution as fine as 1 km. The soil physical texture map was retrieved from a 1:14,000,000 scale map (Institute of Soil Science and Academia Sinica, 1989). The land-use data from the late 1990s were obtained from the land-use datasets (http://www.resdc.cn/), that were developed by the Chinese Academy of Sciences from Landsat Thematic Mapper (TM) scenes, with a spatial resolution of 30 × 30 m. These data were re-sampled to 1 × 1 km using the nearest neighbor resampling method. 3.2.2. Meteorological data With the inverse distance square method, the daily meteorological data (average, maximum and minimum air temperatures, atmospheric pressure, humidity, wind speed, precipitation, and sunshine duration) were interpolated to the entire region with the meteorological data (1999–2010) from 80 meteorological stations in and around the NCP (Fig. 2). 3.2.3. Agricultural data The agricultural data include the measured evapotranspiration and yield data. The evapotranspiration of wheat was measured

by a large-scale weighting lysimeter (2001–2008) from the Luancheng Agricultural Experimental Station (37◦ 55 N, 114◦ 39 E), Yucheng Agricultural Experimental Station (37◦ 00 N, 116◦ 30 E), and Fengqiu Agricultural Experimental Station (35◦ 00 N, 114◦ 30 E) (Fig. 2); these data were obtained from the Chinese Ecosystem Research Network (http://www.cerndata.ac.cn/). The yield data include census data of yield and measured yield. The census data of yield at the county scale of Beijing (2001–2005), Tianjin (2001–2004), Shandong (2001, 2006–2008), Henan (2001–2004), and Hebei (2001–2008) Provinces were derived from the China agricultural annals of 2002–2009 and the provincial agricultural annals. The measured yields of 20 agricultural meteorological stations (Fig. 2) from 2001 to 2010 in the NCP were downloaded from CMDSSS (China Meteorological Data Sharing Service System, http://cdc.cma.gov.cn/index.jsp).

3.2.4. Remote sensing (RS) data The RS data used in this study were during the wheat growing season, including 16-day composite 1-km NDVI from 2000 to 2010, 8-day composite 1 km fPAR from 2000 to 2010, 8-day composite 1-km LAI from 2000 to 2010, and daily land surface temperature from 2000 to 2010 at 1-km resolution. The NDVI data and land surface temperature data were derived from the MOD13A2 and MOD11A1 NASA science data products (http://modis.gsfc.nasa.gov/data/atbd/atbd mod13.pdf), respectively. The fPAR and LAI data (collection 5) were sourced from the MOD15A2 NASA science data product. All of the RS data were distributed by the Land Processes Distributed Active Archive Center (LP DAAC), which is located at the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center (http://lpdaac.usgs.gov). The land cover mask derived from Landsat TM images (http://www.resdc.cn/) was used to eliminate the influence of non-agricultural regions on the RS signal. Consequently, all areas with non-agricultural land were masked out, and RS values were extracted for only those areas containing crops. We also excluded all locations in which the absolute range of NDVI values (maximum–minimum of the annual NDVI) was less than 0.1.

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

The majority of these pixels indicate developed areas. The locations in which the NDVI values were beyond the scope of 0–1 were also excluded because the majority of these pixels indicate surface water areas. Due to cloud contamination, spatially and temporally complete MODIS-fPAR , MODIS-LST and MODIS-LAI were not available during the entire growth stage. The retrieval index, which indicates the number of good retrievals relative to the total number of attempted pixels (Wang et al., 2001), was very low in the study area during the growing season due to the presence of clouds. Non-vegetation and fill value pixels were all excluded based on the QC (quality control) layer for MODIS-fPAR , MODIS-LST and MODIS-LAI. fPAR values derived from both the main radiative transfer method and empirical methods were used, excluding the pixels with geometry problems. LST values and LAI values derived from the main radiative transfer method were used to compare to the simulated results of the model. In addition, the Savitzky–Golay (SG) filter (Chen et al., 2004) was applied to smooth the annual MODIS-fPAR profile and construct high-quality NDVI. All of the geographic information and remote sensing images were processed for projection transformation (Lambert projection). 4. Retrieval implementation As the primary crop in the NCP, winter wheat, which has an area of 30 × l04 km2 , was used as the study subject. We retrieved the Vcmax pattern of winter wheat at a 1-km resolution using the method described in Section 2.2.

325

Fig. 3. The relationship between the Vcmax value retrieved using the assimilation strategy shown in Fig. 1 for 5% randomly selected pixels in the NCP and the corresponding average NDVI value from 2000 to 2010 on the 129th day of the year.

selected pixels was obtained using the global section search algorithm. The results show that all of the functions RMSE = f(Vcmax ) are unimodal, which indicates that the golden section searching algorithm can be used in the optimization design. Based on the MODIS-fPAR from 2000 to 2010, the Vcmax values of sample set A were retrieved using the assimilation strategy in Fig. 1; its relationship with NDVI is shown in Fig. 3.

4.1. Retrieval of the sample set photosynthetic capability 4.2. Interval estimation The pixels in the study area were used as the population, and 5% were randomly selected as the sample (sample set A). The frequency distribution of the average NDVI from 2001 to 2010 on the 129th day (the flowering stage of wheat) of the year showed that the population and sample displayed similar frequency distributions; clearly, the sample had good spatial representation (data not shown). The observed data time window should be scientifically determined as the observed data sequence (MODIS-fPAR ) used to construct the cost function can be the any year from 2001 to 2010. The fPAR sequence from ten years was divided into a calibration period (2001–2007) and a validation period (2008–2010). To determine the best window length, we randomly selected m (m ∈ [1, 7]) year from the calibration period for each simulation. For each set of calibration data, the optimal Vcmax value was determined using the processes in Fig. 1. The optimal Vcmax values were then used in the VIP model to simulate the fPAR , which corresponded to the validation data, and the RMSE between simulation results and corresponding validation data MODIS-fPAR were calculated. For any given m, the average RMSE of each simulation was used for the final evaluation. The above simulation process was performed for the forty randomly selected pixels in the area; the results are shown that a progressive decrease in RMSE was observed between the simulation value and test data as m increased, which indicates that longer time series data observations correlated with improved results. Therefore, we used eleven consecutive years (2000–2010) of MODIS-fPAR data for parameter retrieval in this study. The golden section search algorithm used in this study is a one-dimensional search algorithm with the precondition that the optimization function is unimodal. Because the ecosystem model was used as the optimization function in this study, testing for unimodality is difficult; therefore, the Vcmax value of the 100 randomly

The Vcmax values in the similar environment (which had a similar NDVI value) were located in a narrower interval than the interval determined by its special attribution, as shown in Fig. 3. For example, in the drought condition, NDVI is usually lower than 0.4, and its corresponding Vcmax value is located in the interval (40 ␮mol C m−2 s−1 , 60 ␮mol C m−2 s−1 ) which is narrower than the interval determined by its special attribution (40 ␮mol C m−2 s−1 , 120 ␮mol C m−2 s−1 ). Therefore, sample set A was divided into 50 sub-sample sets {A1 , A2 , A3 , . . ., A50 } based on the NDVI value to narrow the range of variation in the parameter. Although the number of sub-sample sets varied, they numbered greater than 30, which indicates that {A1 , A2 , A3 , . . ., A50 } are statistically large samples. The interval estimation method was used for each sub-sample set An (n ∈ [1, 2, 3, . . ., 50]) to construct the population’s confidence interval of Vcmax . The precondition for applying the interval estimation is a random sample with a normal distribution. Therefore, the Kolmogorov–Smirnov test (K–S test) (Justel et al., 1997) was performed to determine whether the sub-sample set An (n ∈ [1, 2, 3, . . ., 50]) from the population had a normal distribution. The results show that sub-sample An was from a population with a normal distribution when its corresponding NDVI is greater than 0.3. Therefore, the interval estimation method can be used for the subsample set An (n ∈ [15, 16, 17, . . ., 50]), with a corresponding NDVI greater than 0.3. When the pixel NDVI is less than 0.3, the photosynthetic parameter does not obey a normal distribution, and the photosynthetic parameter is less than 50 ␮mol C m−2 s−1 (Fig. 3); therefore, the new search interval is defined as [40 ␮mol C m−2 s−1 , 50 ␮mol C m−2 s−1 ] in this study. Based on Eqs. (7) and (8), the population mean’s 95% confidence interval and the population standard deviation’s 95% confidence interval were calculated using each sub-sample set An (n ∈ [15, 16,

326

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

Fig. 4. Interval estimation for the (a) population mean and (b) population SD; the filled dots and empty circles connected by a solid line are the upper and lower confidence limits for the 95% confidence interval of the population mean (in panel a)/population SD (in panel b) for each sub-sample set, the solid line and the dash line are the fit lines for the upper and lower confidence limits for the 95% confidence interval, respectively.

17, . . ., 50]); both the upper and lower limits of the population mean’s 95% confidence interval and the population standard deviation’s 95% confidence interval have good relationship with the NDVI, as shown in Fig. 4. Eq. (9) were fit to the relationship between the population mean’s confidence interval and its corresponding NDVI, Eq. (10) were fit to the relationship between the population standard deviation’s confidence interval and its corresponding NDVI.

⎧ ⎪ ⎨ Vc ⎪ ⎩ Vc 

up

= 40 +

low

= 40 +

72 1 + 10(3.96−8.5x) 67

(9)

1 + 10(4.62−9.7x)

up = −0.13 + 13.8 × low = −0.08 + 9.1 ×



2



2

/2 × e−3.47(x−0.48)

Taihang piedmont, Henan and Jiangsu, which accounted for 10% of the pixels. The pixels with photosynthetic capabilities lower than 50 ␮mol C m−2 s−1 were mostly in the eastern coastal zone, which has alkila-saline soil, and the southeast portion of the NCP; these areas accounted for 25% of the pixels. 5.2. Validation using the crop yield The census data on county-level grain yield show wide range variability that coincides with the heterogeneity of the soil, water resources and climatic conditions at the study site. The spatial variation of the average simulated yield from 2001 to 2008 (Fig. 6a) was highly consistent with the yield census data; the R2 was 0.63,

(10)

/2 × e−3.47(x−0.48)

where Vc up and Vc low are the upper and lower limits of the population mean’s 95% confidence interval,  up and  low are the upper and lower limits of the 95% confidence interval for the population standard deviation, x is the average NDVI from 2000 to 2010 on the 129th day of the year. The population parameter’s standard deviation ranged from 0 to 12 ␮mol C m−2 s−1 . In the area with a relatively stable environment (primarily in areas with rich soil and sufficient irrigation or areas with heavy alkila-saline soil and insufficient irrigation), the standard deviation was low. This result indicates that the data points tend to be close to the mean, and the Vcmax value has a narrow variation range. High standard deviation values were typically observed in areas with large environmental changes, where the NDVI values were between 0.45 and 0.6. 5. Method validation 5.1. Spatial variation of the photosynthetic pattern Based on the new parameter searching interval ([Vc low − 1.96 up , Vc up + 1.96 up ]), the NCP winter wheat Vcmax pattern was generated using the golden section search algorithm, as shown in Fig. 5. The Vcmax value ranged from 40 to 120 ␮mol C m−2 s−1 with an average value 70 ␮mol C m−2 s−1 and a SD 16 ␮mol C m−2 s−1 ; it decreased in a south-to-north direction. The pixels with photosynthetic capabilities greater than 90 ␮mol C m−2 s−1 were mostly in

Fig. 5. Spatial pattern of Vcmax with 1-km resolution in the NCP (the blank area is not farmland).

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

327

Fig. 6. Spatial variation of (a) the average simulated winter wheat yield from 2000 to 2010 and (b) average cumulative ET during wheat growing seasons from 2000 to 2010.

and the RMSE% was 18.3% (the RMSE% equals the RMSE divided by the observed data). The average yield in the study area was 0.59 kg m−2 (with a 0.09 kg m−2 SD), which was similar to the average observed yield (0.64 kg m−2 with the SD 0.086 kg m−2 ) of the 149 samples in the study area (Liu et al., 2011). The R2 of the simulated yield and census yield data from 2001 to 2008 (eight datasets for 8 years) ranged from 0.56 to 0.70, and the RMSE ranged from 0.08 to 0.13 kg m−2 (the RMSE% ranged from 13 to 25%) (Fig. 7a). The results indicate that the model can perform satisfactorily in a yield pattern simulation at the county level. The validation results between the simulated and measured yield at 20 agricultural meteorological stations from 2001 to 2010 show that the R2 between the simulated and measured yield range from 0.53 to 0.97, and the RMSE ranges from 0.034 kg m−2 to

0.15 kg m−2 (the RMSE% ranges from 6.1 to 25%). The above-ground biomass in 2011 and 2007 was validated using measured field data collected in Tongzhou and Luancheng, respectively. The simulated above-ground biomass was consistent with the measured values in Tongzhou (R2 = 0.93) and Luancheng (R2 = 0.97). The RMSE% between the simulated and measured above-ground biomass in Tongzhou (28%) was higher than in Luancheng (13%).

5.3. Validation using the leaf area index The measured LAI at the Tongzhou and Luancheng agricultural stations and MODIS-LAI in the NCP were both used for model validation. The simulated LAI was consistent with the measured values in Tongzhou and Luancheng (Fig. 8). The RMSE% between the

Fig. 7. Comparison of simulated yield using the retrieved Vcmax value with (a) census data on the yield at the county level from 2001 to 2008 and (b) the measured yield at 20 agricultural meteorological stations from 2001 to 2010.

328

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

Fig. 8. Comparison of simulated LAI with measured field data in (a) the Luancheng agricultural station in 2007 and (b) the Tongzhou agricultural station in 2011.

Table 1 R2 and RMSE between the simulated LAI and MODIS-LAI.

5.4. Validation using evapotranspiration data

Year

R2

RMSE

n

Year

R2

RMSE

n

2001 2002 2003 2004 2005

0.19 0.31 0.38 0.26 0.24

2.18 2.08 2.09 2.05 1.67

325,788 317,206 324,900 325,556 325,633

2006 2007 2008 2009 2010

0.28 0.24 0.26 0.28 0.44

2.02 2.18 1.87 2.03 2.11

325,591 325,715 325,352 325,450 325,469

n is the number of valid MODIS-LAI data retrieved from the main algorithm (QC < 64).

simulated and measured LAI in Tongzhou (19%) was lower than in Luancheng (14%). The valid MODIS-LAI data from 2001 to 2010 retrieved from the main algorithm (QC < 64) were used in the LAI validation. Because the 129th day was during the flowering stage of wheat, and the LAI value reach a maximum near the 129th day, the simulated LAI values for the 129th day in 2001 to 2010 were compared with the corresponding MODIS-LAI values (Table 1). The R2 and RMSE between the simulated and measured LAI ranged from 0.19 to 0.44 and 1.87 to 2.11, respectively.

The spatial variation of the average simulated ET from 2000 to 2010 is shown in Fig. 6b. The cumulative ET during wheat growing seasons ranged from 280 to 445 mm with the average value 380 mm and 28 mm SD; it decreased in a south-to-north direction. The simulated monthly and accumulated ET values for 10 days from 2001 to 2008 were compared with the measured data from Luancheng, Fengqiu and Yucheng (Fig. 9). The ET simulated on monthly and 10-day scales were both consistent with the field data and had R2 values greater than 0.79 (monthly scale) and 0.72 (10-day scale) as well as RMSE values less than 18.2 mm (monthly scale) and 4.9 mm (10-day scale) for the three sites.

5.5. Validation using the land surface temperature The recently released MODIS LST V5 products (Wan, 2008; Coll et al., 2009), which provide global coverage and high resolution (the finest is a 1-km grid), are a powerful tool for improving our regionalscale energy budget study. The spatial variation of the time-series correlation coefficient and RMSE calculated for the simulated

Fig. 9. Comparison of the simulated cumulative ET and observed field data using a lysimeter at the Luancheng, Fengqiu and Yucheng agricultural stations from 2001 to 2008 (a represents the cumulative ET on a 10-day scale, and b represents the cumulative ET on a monthly scale).

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

329

Fig. 10. The spatial variation of a time-series (a) correlation coefficient and (b) RMSE calculated for the simulated LST and MODIS-LST at a 1-km resolution from 2001 to 2010.

LST and MODIS-LST from 2001 to 2010 is shown in Fig. 10. The average values of the correlation coefficient and RMSE for all of the pixels in the NCP were 0.92 and 2.58 ◦ C, with SDs of 0.02 and 0.13 ◦ C, respectively. Correlation coefficients in greater than 90% of the pixels were higher than 0.9 with the maximum frequency 0.93, and the RMSE values for greater than 90% of the pixels were less than 5.0 ◦ C with the maximum frequency at 2.8 ◦ C, which indicates that the model performs well when it is used to simulate the LST time-variant process. Due to cloud cover, seasonal snow and instrument problems, the current LST results are spatial and temporally discontinuous.

Therefore, the LST simulated on the 129th day of 2008, which provided more data from the main radiative transfer method, was used for the comparison. A considerably similar spatial pattern was observed for the simulated LST and MODIS-LST; the land surface temperature decreased from west to east (Fig. 11). The simulated LST ranged from 16.5 to 33.0 ◦ C with an average value and standard deviation of 25.8 ◦ C and 1.3 ◦ C, respectively. The MODIS-LST ranged from 17.0 to 39.5 ◦ C with an average value and standard deviation of 27.4 ◦ C and 2.7 ◦ C, respectively. The simulated LST for 90% of the pixels was lower than 28 ◦ C, while the MODIS-LST for 60% of the pixels was lower than 28 ◦ C.

Fig. 11. Spatial variation of (a) the simulated LST and (b) MODIS-LST on the 129th day of 2008.

330

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

6. Discussion 6.1. Factors impact spatial pattern of Vcmax Rubisco synthesis and the Vcmax value are closely related to nitrogen concentration (Bruck and Guo, 2006; Imai et al., 2008) and water stress (Flexas et al., 2006). The area with sufficient irrigation and fertile soil in the north of NCP has photosynthetic capabilities greater than 90 ␮mol C m−2 s−1 . In contrast, poor soil conditions and insufficient irrigation may limit crop photosynthetic capabilities. The area with insufficient irrigation and alkila-saline soil are mostly in the southeast portion of the NCP, and has photosynthetic capabilities lower than 50 ␮mol C m−2 s−1 . This is similar to reports by Alonso et al. (2009), wherein the Vcmax of wheat was 93.2 ␮mol C m−2 s−1 in areas with sufficient irrigation and nutrition, and Kothavala et al. (2005), wherein the photosynthetic value for wheat grown in silt–clay soil with the daily temperature range −3.7 to 33.9 ◦ C and the annual precipitation 835 mm was 93 ␮mol C m−2 s−1 . Certain mechanical investigations have noted that nitrogen use efficiency in non-saline soil exceeds efficiency in saline soil by 15% (Irshad et al., 2008), and N uptake is independent of the N rates and is mainly influenced by soil salinity at high salinities (Chen et al., 2010). These reports demonstrate that soil alkalization reduces nitrogen availability (Li et al., 2008a) due to salt osmosis. If the leaf nitrogen levels decreased to 1.19 g m−2 , the Vcmax would correspondingly decrease to 58 ␮mol C m−2 s−1 (Muller et al., 2005). Consequently, the average photosynthetic values in the inland alkila-saline soil (54 ␮mol C m−2 s−1 ) and costal alkilasaline soil (44 ␮mol C m−2 s−1 ) were lower than in non-alkila-saline soil (68 ␮mol C m−2 s−1 ), sandy soil (64 ␮mol C m−2 s−1 ) and sajong black soil (78 ␮mol C m−2 s−1 ). 6.2. Errors in yield, LAI, ET and LST simulation Errors in yield simulation are mainly attributed to irrigation amount and planting area. The irrigation level was unavailable in the study area and set at 240 mm per year. This value may be underestimated in certain highly productive plots and overestimated in certain low productivity plots, which could lead to an underestimation or overestimation in these plots, respectively, as shown in Fig. 7b. The bias in planting area estimates may result in uncertainty for yield. Crop yield is better modeled in pixels where the wheat field percentage is greater. For example, in the agricultural meteorological station Shangqiu (Fig. 2a), which has a wheat field percentage greater than 98%, the R2 between the simulated and measured yield is 0.97, and the RMSE is 0.054 kg m−2 (the RMSE% is 9.1%); while in agricultural meteorological station Huangye (Fig. 2a), which has a wheat field percentage less than 30%, the R2 between the simulated and measured yield is 0.53, and the RMSE is 0.071 kg m−2 (the RMSE% is 18%). Compared with yield measured in 236 field plots in the NCP, Xiong (2009) showed an R2 of 0.52, and the RMSE% was 16.6% between simulated wheat yield using a CERES-wheat model and measured data. Although the uncertainties for the current MODIS-LAI products (±1.0) do not satisfy the GCOS accuracy requirement (±0.5), they provide an acceptable description of the global LAI spatial-temporal variation (Fang et al., 2012). The R2 between the simulated and measured LAI ranged from 0.19 to 0.44. By comparing the global MODIS-LAI with the measured data, Fang et al. (2012) found that the R2 and RMSE values between the measured LAI and MODISLAI were 0.165 and 1.24, respectively, which was similar to our reports herein. Considering the MODIS-LAI underestimation, the LAI simulation results are acceptable. Although the simulated ET was quite close to the observation data, the deviations still exist in some period. They may be related

to the errors of the model and the observation. Sensitivity analysis illustrates that ET is quite sensitive to soil resistance parameters contributing to the soil evaporation, the uncertainties of annual ET were 16.5% when the key parameter are randomly sampled in their uncertainty ranges (Mo et al., 2012). Except for the error in model parameters, there are also errors in observation. Measurements from the large weighing lysimeter could be inflated due to the “oasis effect” (Jiang et al., 2008), while the eddy covariance technique does not always achieve energy budget closure, and this may introduce the underestimation of latent and sensible heat flux (Liu et al., 2009). Therefore, the simulated ET was slightly lower than the corresponding observed lysimeter data at the three sites (the slop of the regression is less than 1 in Fig. 9) and the corresponding observed lysimeter data from other reports. For example, the average ET measured by lysimeter in the wheat growing seasons was 453 mm in Luancheng from 1995 to 2000 (Liu et al., 2002) and 458 mm in Yucheng from 1986 to 2007 (Liu et al., 2009), while it was 436 mm in Luancheng and 440 mm in Yucheng in this study. Compared with the observed data from the eddy covariance technique, the simulated result had a higher cumulative ET. For example, the simulated total ET during the wheat growing seasons in Yucheng in 2006 was 435 mm, while it was 412 mm as determined using the eddy covariance technique (Lei and Yang, 2010). Mo et al. (2005) reported that the ET during the wheat growing season ranged from 330 to 500 mm in the NCP. In this study, the simulated ET and observed data in three sites from 2001 to 2008 ranged from 410 to 480 mm during the wheat growing seasons, which demonstrated that the observation data and the simulation results were both in an acceptable range. Comparisons between V5 LSTs and in situ values in 47 clear-sky cases indicate that the accuracy for the MODIS-LST products was better than 1 K in most cases (39 out of 47), and the RMSE was less than 0.7 K for all 47 cases (Wan, 2008). The RMSE between the MODIS-LST and simulated LST in this study was similar to reports by Wang et al. (2009) in the upper Tone River Basin (2.22 K). These results suggest that the model can simulate regional LST patterns well; however, it underestimated the LST in certain high-altitude areas. Underestimates of the modeled daytime LSTs were generally restricted to areas with the highest altitudes; these errors are potentially attributed to errors in the input air temperature. Because no meteorological sites were available, a homogeneous lapse rate (6.5 K km−1 ) was used to account for the air temperature dependence on elevation; therefore, errors in the interpolated air temperature may have occurred in the uppermost mountain regions (Wang et al., 2009). Meanwhile, irrigation can reduce the land surface temperature; however, because the exact irrigation levels were unavailable in the study area, the average irrigation level 240 mm per year was used. This value may be slightly higher than the actual levels in the area at high altitudes, which may have also contributed to the LST underestimation. 6.3. Comparison with other Vcmax retrieval approaches Through calibrating the parameters, initializing model variables/parameters, adjusting the model state variables with the retrieved remote sensing variables, assimilating observed remotely sensed data improved the ecosystem model predictions at the regional scale (Dorigo et al., 2007). Compared with replacing model state variables with RS data, using the RS data adjust model parameter is more flexible for remote sensing data assimilation because its use is unrestricted by using real-time remote sensing data. Many approaches have been described for retrieving biophysical parameters from remote sensing data. For Vcmax , the regional winter wheat Vcmax pattern in the NCP was derived from the model’s predicted photosynthetic capacity through regression of the accumulated NDVI during the growing season and grain yield at the

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

331

accumulated NDVI during the growing season and grain yield at the county level (it is similar to a statistical approach); this practice may misestimate the values of certain pixels due to the NDVI fusion. However, the problem with fusion of the NDVI at the county level can be avoided to a certain extent when each pixel is considered in determining Vt . Although the method utilized in the current study was complex when compared with the regression method proposed by Hu and Mo (2012), it only required remote sensing data to obtain the Vcmax pattern, whereas remote sensing data and census data were both used to retrieve regional Vcmax pattern in the regression method. Considering the difference observed between the two methods were very small, the method utilized in the current study will show more advance in the census data sparse region. 6.4. Vcmax issues

Fig. 12. The relationship between the Vcmax value retrieved by the two methods and its corresponding average NDVI value from 2001 to 2010 on the 129th day of the year, Vt is the Vcmax value obtained using the method described in the current study, while Vj is the Vcmax value obtained using the method introduced by Hu and Mo (2012).

county level (Hu and Mo, 2012, details can be seen in Appendix C). Although the Vcmax spatial variation obtained in this study was highly consistent with Hu and Mo (2012), certain differences were observed between the two patterns. The Vcmax values for 1500 pixels in the study area were obtained using the method described in this study (referred to as Vt ) as well as the method introduced by Hu and Mo (2012) (referred to as Vj ) and were randomly selected and compared to illustrate the differences. Vt was highly related to Vj , with R2 of 0.85. The differences between Vt and Vj (namely Vt − Vj ) were approximately ±30 ␮mol C m−2 s−1 ; 90% of the pixels had differences less than or equal to ±15 ␮mol C m−2 s−1 (Fig. 12). When the pixel NDVI was between 0.3 and 0.45, which indicates that the pixel had a relatively low photosynthetic capability, Vj was higher than Vt , whereas the opposite trend was observed in pixels with relatively high photosynthetic capabilities (the NDVI was between 0.5 and 0.65). Compared with the measured data, the R2 between the measured above-ground biomass, measured LAI and corresponding simulated data with Vt (such as for Luancheng) was slightly higher than with Vj (such as for Tongzhou), but the RMSE between the measured and simulated data with Vt was slightly lower than with Vj . This result was observed because Vj was obtained by regressing the

The Vcmax value can be affected by any factors that affect the level of activated enzyme (Rubisco). Its average value is affected by the extent of environmental stress, such as available nitrogen (Makino, 2003; Bruck and Guo, 2006), heavily salinized soil and groundwater (Tester and Davenport, 2003; James et al., 2006). These environmental stresses will result in inter-annual Vcmax pattern variations. Hence, using a fixed Vcmax pattern may produce errors in the yield and ET simulations (Hu and Mo, 2014). One hundred twenty pixels were randomly selected in the study. For each pixel, the yearly Vcmax values (referred to as unfixed Vcmax ) from 2001 to 2010 were retrieved using the MODIS-fPAR for the corresponding year with the method utilized in the current study. Therefore, the Vcmax values changed from year to year, which indicates the inter-annual variation in environmental influences. The simulated yield and ET from 2001 to 2010 with unfixed Vcmax are referred to as Yo and ETo , and those with fixed Vcmax , as determined using the MODIS-fPAR from 2001 to 2010 in the current study, are referred to as Ym and ETm . The difference between Yo and Ym as well as between ETo and ETm decreased with the enhanced crop photosynthetic capability (Fig. 13). Usually rain-fed fields and saline–alkaline soil with relatively low Vcmax value has larger environment changes than irrigated fields with relatively high Vcmax value (Hu and Mo, 2014), using a fixed Vcmax pattern is lack of concern for the expression of the inter-annual environmental change, therefore the difference between Yo and Ym as well as between ETo and ETm decreased with the increase of Vcmax value. The differences between the Yo and Ym as well as between ETo and ETm include differences due to the Vcmax difference and sampling error. The yield and ET differences induced by the Vcmax

Fig. 13. The RMSE between (a) simulated yield and (b) simulated ET with a fixed Vcmax and an unfixed Vcmax for different Vcmax levels, the squares in the boxes correspond to the mean value from 2001 to 2010, and the whiskers extend to the minimum and maximum values. A, B, C, and D are the pixels with Vcmax values ranging from 40 to 60 ␮mol C m−2 s−1 , 60 to 80 ␮mol C m−2 s−1 , 80 to 100 ␮mol C m−2 s−1 and 100 to 120 ␮mol C m−2 s−1 , respectively.

332

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

change cannot be separated through the direct comparisons above. With respect to the large spatial variation in climate and soil in the NCP, the t-test method (DeGroot and Schervish, 2012) was used to test for significance differences between Yo and Ym as well as between ETo and ETm . The t-values between Yo and Ym (ranged from 0.391 to 1.975) as well as between ETo and ETm (ranged from 0.889 to 1.957) from 2001 to 2010 at four Vcmax value levels were all lower than the t0.05 threshold (2.021), which indicates no significant differences between Yo and Ym or between ETo and ETm . In other words, the differences between Yo and Ym as well as between ETo and ETm were mainly due to sampling error (P < 0.05). Although the Vcmax pattern depends on the near real-time remote sensing data, using a fixed Vcmax pattern to predict yield and ET over ten years in the NCP is acceptable based on the aforementioned significance tests.

where Rn , H, and LE are the absorbed net radiation (W m−2 ), sensible heat flux (W m−2 ) and latent heat flux (W m−2 ), respectively. The subscripts “c” and “g” refer to the canopy and the soil surface, respectively. G is the soil heat flux (W m−2 ). Cc and Cg are the effective heat capacity (J m−2 K−1 ) of the canopy and soil, respectively. The equations for calculating Rn , H, LE, Cc and Cg can be found in Mo and Liu (2001) and Mo et al. (2012). In the MODIS LST product, LST is the radiometric (kinetic) temperature related to the thermal infrared (TIR) radiation emitted from the land surface observed by an instantaneous MODIS observation (Wan, 2008), where the land surface is the top of the canopy in vegetated areas or the soil surface in bare areas. For a model grid of mixed vegetation and bare soil, LST can be estimated from Tc and Tg if the emissivity of the model grid is assumed as homogeneous (Becker and Li, 1995; Norman and Becker, 1995)

7. Conclusions Tsim = The approach developed in this paper introduces the interval estimation method to a one-dimensional searching algorithm and offers a solution for assimilating remote sensing data with the VIP ecosystem model at regional scale. The pattern of photosynthetic parameter – Vcmax for winter wheat in the NCP was obtained with the developed approach. Satisfactory agreement was obtained between the simulated and measured data with respect to yield, evapotranspiration, surface temperature and leaf area index, which demonstrates that the Vcmax pattern is reliable. The approach developed in the current study is unconstrained by grain census data; it shows more advance in the census data sparse region when compared with a statistical photosynthetic parameter retrieval approach. However, the Vcmax pattern retrieved by the approach depends on near real-time remote sensing data and is a time-invariant pattern for ET and yield simulation over ten years. Considering the simulated yield and ET with a time-invariant Vcmax value showed no significant differences from time-variant Vcmax values, using a time-invariant Vcmax pattern in 10-year projections for yield and ET in the NCP is acceptable. Acknowledgements This study was jointly supported by the Natural Science Foundation of China (NSFC) grants (31171451, 31300374), the Chinese Ministry of Science and Technology Projects (2010CB428404), the Key Project for the Strategic Science Plan in IGSNRR, CAS (2012ZD003), the Natural Science Foundation of China (41171448) and the National Science Foundation for Post-doctoral Scientists of China (07Z76034Z1). We thank the anonymous reviewers and the editors of the journal for their comments that improved an earlier draft of this paper. Appendix A. Land surface energy balance Energy fluxes are described with a two-source scheme that individually discerns the canopy and the soil surface. The transfer of solar radiation within the crop canopy is simulated with 20 sublayers, distinguishing the direct and diffuse components of the visible (300–700 nm) and near-infrared (700–1300 nm) fractions (Mo and Liu, 2001). The temperatures of the canopy (Tg ) and soil surface (Tc ) are governed by the two-sources soil-canopy energy balance sub-model (Sellers et al., 1996): Cc

∂Tc = Rnc − Hc − LEc ∂t

(A1)

Cg

∂Tg = Rng − Hg − LEg − G ∂t

(A2)

[(1 − )εc Tc4 + εg Tg4 ]

1/4

ε˛

(A3)

where ε˛ is the apparent emissivity of the canopy (taken as 0.98), εc and εg are the emissivities of the soil and canopy, taken as 0.96 and 0.98, respectively.  is the sky view factor for soil and can be expressed as follows for most crops:



 = exp −0.8





LAI

(A4)

Appendix B. Crop dynamics The daily growth of the crop dry mass, denoted as M (mol C m−2 d−1 ), is expressed as the balance of gross photosynthesis, respiration and senescence in the vegetative stage, namely: dMi = ˛i Ac − Rci − Rmi − Di dt

(B1)

where the subscript i represents either the stem (s), leaf (l) and root (r), Ac is the daily net photosynthesis (mol C m−2 d−1 ) at the canopy level determined by integrating the instantaneous An values over daytime, ˛ is the partition coefficient, Rc (mol C m−2 d−1 ) is the growth respiration, Rm (mol C m−2 d−1 ) is the maintenance respiration, and D (mol C m−2 d−1 ) is the senescence rate. The daily growth of the LAI is determined by the daily dry mass allocated to the leaf and the specific leaf area (SLA m2 g c−1 ) (Eq. B2). Based on the results from the field tests in the study area, the specific leaf area is set as 0.07. LAI = SLA × Ml

(B2)

Appendix C. Vcmax retrieval method proposed by Hu and Mo (2012) (1) Seeking the relationship between NDVI and crop yield The linear relationship between yield and the NDVI is established by the cumulative NDVI data and the census data of yield. The census data is the county-level wheat yield data of the 104 counties in the study region in 2001, and the cumulative NDVI data is the corresponding county-level average value of cumulative NDVI from March to June 2001. The regressive equation is expressed as follows: yyield = 747.71xNDVI + 2371.40

(C1)

where yyield is the county-level yield of winter wheat, and xNDVI is the county-level average value of accumulative NDVI. (2) Seeking the relationship between Vcmax and the crop yield by the VIP model Different Vcmax is used to drive the VIP model to obtain the simulated yield in Shijiazhuang, which is a typical area with wheat-maize rotation system in the NCP. The Vcmax and the

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

corresponding simulated yield data are used to establish the relationship between Vcmax and yield, the regressive equation is expressed as follows: yVc max = 17.46e0.0003xyield

(C2)

where yVc max is the value of Vcmax , and xyield is the winter wheat yield. (3) Obtaining the relationship between Vcmax and NDVI with the aid of crop yield. Since the relationship between Vcmax and yield is controlled by the mechanism of crop growth, the relationship between the yield and Vcmax in Shijiazhuang can be scaled up to the study region. Based on the assumption that the VIP model has good performance in yield simulation at a regional scale, which means that the simulated yield is equal to the census data of yield, the relationship between Vcmax and NDVI can be derived from the two equations above and is expressed as follows: yVc max = 17.46e(0.23xNDVI +0.71)

(C3)

where yVc max is the value of Vcmax at each grid and xNDVI is the value of the NDVI at each grid. References Alonso, A., Perez, P.A., Martinez-Carrasco, R., 2009. Growth in elevated CO2 enhances temperature response of photosynthesis in wheat. Physiol. Plant. 135, 109–120. Atzberger, C., Jarmer, T., Schlerf, M., Koetz, B., Werner, W., 2003. Retrieval of wheat bio-physical attributes from hyperspectral data and SAILH + PROSPECT radiative transfer model. In: Proceedings of the Third EARSeL Workshop on Imaging Spectroscopy, Herrsching, Germany. Bastiaanssen, W.G.M., Ali, S., 2003. A new crop yield forecasting model based on satellite measurements applied across the Indus Basin, Pakistan. Agric. Ecosyst. Environ. 94 (3), 321–340. Basnyat, P., Mcconkey, B., Lafond, G.R., Moulin, A., Pelcat, Y., 2004. Optimal time for remote sensing to relate to crop grain yield on the Canadian prairies. Can. J. Plant Sci. 84, 97–103. Becker, F., Li, Z., 1995. Surface temperature and emissivity at various scales: definition, measurement and related problems. Remote Sens. Rev. 12, 225–253. Becker-Reshef, I., Vermote, E., Lindeman, M., Justice, C., 2011. A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data. Remote Sens. Environ. 114, 1312–1323. Bruck, H., Guo, S.W., 2006. Influence of N form on growth and photosynthesis of Phaseolus vulgaris L. plants. J. Plant Nutr. Soil Sci. 169, 849–856. Centritto, M., Loreto, F., Chartzoulakis, K., 2003. The use of low CO2 to estimate diffusional and non-diffusional limitations of photosynthetic capacity of saltstressed olive sap lings. Plant Cell Environ. 26, 585–594. Chen, J., Jonsson, P., Tamura, M., Gu, Z., Matsushita, B., Eklundh, L., 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. Remote Sens. Environ. 913, 332–344. Chen, W.P., Hou, Z.A., Wu, L.S., Liang, Y.C., Wei, C.Z., 2010. Effects of salinity and nitrogen on cotton growth in arid environment. Plant Soil. 326, 61–73. Choi, B.Y., Roh, K.S., 2003. UV-B radiation affects chlorophyll and activation of Rubisco by Rubisco activase in Canavalia ensiformis L. leaves. J. Plant Biol. 46 (2), 117–121. Coll, C., Hook, S.J., Galve, J.M., 2009. Land surface temperature from the advanced along-track scanning radiometer: validation over inland waters and vegetated surfaces. IEEE Trans. Geosci. Remote Sens. 47, 350–360. Dawson, T.P., North, P.R.J., Plummer, S.E., Curran, P.J., 2003. Forest ecosystem chlorophyll content: implications for remotely sensed estimates of net primary productivity. Int. J. Remote Sens. 24 (3), 611–617. Dente, L., Satalino, G., Mattia, F., Rinaldi, M., 2008. Assimilation of leaf area index derived from ASAR and MERIS data into CERES-wheat model to map wheat yield. Remote Sens. Environ. 112, 1395–1407. DeGroot, M.H., Schervish, M.J., 2012. Probability and Statistics, fourth edition. Pearson Education, Inc. Doraiswamy, P.C., Hatfield, J.L., Jackson, T.J., Akhmedov, B., Prueger, J., Stern, A., 2004. Crop condition and yield simulations using Landsat and MODIS. Remote Sens. Environ. 92 (4), 548–559. Dorigo, W.A., Zurita-Milla, R., de Wit, A.J.W., Brazile, J., Singh, R., Schaepman, M.E., 2007. A review on reflective remote sensing and data assimilation techniques for enhanced agroecosystem modeling. Int. J. Appl. Earth Obs. Geoinf. 9, 165–193. Farquhar, G.D., Von, C.C., Berry, J.A., 1980. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149, 78–99. Fang, H., Liang, S., Hoogenboom, G., Teasdale, J., Cavigelli, M., 2008. Corn-yield estimation through assimilation of remotely sensed data into DSSAT-CERES. Int. J. Remote Sens. 29 (10), 3011–3032. Fang, H., Wei, S., Liang, S., 2012. Validation of MODIS and vegetation LAI products using global field measurement data. Remote Sens. Environ. 119, 43–54.

333

Flexas, J., Ribas-Carbo, M., Bota, J., Galmes, J., Henkle, M., artinez-Canellas, S., Medrano, H., 2006. Decreased Rubisco activity during water stress is not induced by decreased relative water content but related to conditions of low stomatal conductance and chloroplast CO2 concentration. New Phytol. 172 (1), 73–82. Gonzalez-Real, M.M., Baille, A., 2000. Changes in leaf photosynthetic parameters with leaf position and nitrogen content within a rose plant canopy (Rosa hybrida). Plant Cell Environ. 23, 351–363. Haboudane, D., Miller, J.R., Tremblay, N., Zarco-Tejada, P.J., Dextraze, L., 2002. Integrated narrow-band vegetation indices for prediction of crop chlorophyll content for application to precision agriculture. Remote Sens. Environ. 81 (2–3), 416–426. Hazarika, M.K., Yasuoka, K., Ito, A., Dye, D., 2005. Estimation of net primary productivity by integrating remote sensing data with an ecosystem model. Remote Sens. Environ. 94, 298–310. Houborg, R., Cescatti, A., Migliavacca, M., Kustas, W.P., 2013. Satellite retrievals of leaf chlorophyll and photosynthetic capability for improved modeling of GPP. Agric. For. Meteorol. 177, 10–23. Hu, S., Mo, X., 2012. Interpreting spatial heterogeneity of crop yield with a process model and remote sensing. Ecol. Model. 222, 2530–2541. Hu, S., Mo, X., 2014. Prediction of crop productivity and evapotranspiration with two photosynthetic parameter regionalization methods. J. Agric. Sci. 152 (1), 119–133. Huang, Z., Turner, B.J., Dury, S.J., Wallis, I.R., Foley, W.J., 2004. Estimating foliage nitrogen concentration from HYMAP data using continuum removal analysis. Remote Sens. Environ. 93 (1–2), 18–29. Huang, L., Ning, Z.Y., Zhang, X.L., 2010. Impact of caterpillar disturbance on forest net primary production estimation in China. Ecol. Indic. 10 (6), 1144–1151. Imai, K., Suzuki, Y., Mae, T., Makino, A., 2008. Changes in the synthesis of Rubisco in rice leaves in relation to senescence and N influx. Ann. Bot. 101, 135–144. Institute of Soil Science, Academia Sinica (compiled), 1989. The Soil Atlas of China. Cartographic Publishing House, Beijing, China. Irshad, M., Eneji, A.E., Yasuda, H., 2008. Comparative effect of nitrogen sources on maize under saline and non-saline conditions. J. Agron. Crop Sci. 194, 256–261. James, R.A., Munns, R., Von Caemmerer, S., Trejo, C., Miller, C., Condon, T., 2006. Photosynthetic capacity is related to the cellular and subcellular partitioning of Na+ , K+ and Cl− in salt-affected barley and durum wheat. Plant Cell Environ. 29, 2185–2197. Jiang, J., Zhang, Y., Wegehenkel, M., Yu, Q., Xia, J., 2008. Estimation of soil water content and evapotranspiration from irrigated cropland on the North China Plain. J. Plant Nutr. Soil Sci. 171, 751–761. Justel, A., Pena, D., Zamar, R., 1997. A multivariate Kolmogorov–Smirnov test for goodness of fit. Stat. Probabil. Lett. 35 (3), 251–259. Kattge, J., Knorr, W., 2007. Temperature acclimation in a biochemical model of photosynthesis: a reanalysis of data from 36 species. Plant Cell Environ. 30 (9), 1176–1190. Kothavala, Z., Arain, M.A., Black, T.A., Verseghy, D., 2005. The simulation of energy, water vapor and carbon dioxide fluxes over common crops by the Canadian Land Surface Scheme (CLASS). Agric. For. Meteorol. 133, 89–108. Launay, M., Guerif, M., 2005. Assimilating remote sensing data into a crop model to improve prediction performance for spatial application. Agric. Ecosyst. Environ. 111 (1–4), 321–339. Lei, H.M., Yang, D.W., Shen, Y.J., Liu, Y., Zhang, Y.C., 2011. Simulation of evapotranspiration and carbon dioxide flux in the wheat-maize rotation croplands of the North China Plain using the Simple Biosphere Model. Hydrol. Process. 25 (20), 3107–3120. Lei, H.M., Yang, D.W., 2010. Interannual and seasonal variability in evapotranspiration and energy partitioning over an irrigated cropland in the North China Plain. Agric. For. Meteorol. 150, 581–589. Leuning, R., 1997. Scaling to a common temperature improves the correlation between the photosynthesis parameters J(max) and V-cmax. J. Exp. Bot. 48 (307), 345–347. Li, W.Q., Li, X., Khan, M.A., Gul, B., 2008a. Relationship between soil characteristics and halophytic vegetation in coastal region of North China. Pak. J. Bot. 40, 1081–1090. Liu, X.Y., He, P., Jin, J.Y., Zhou, W., Sulewski, G., Phillips, S., 2011. Yield gaps, indigenous nutrient supply, and nutrient use efficiency of wheat in China. Agron. J. 103 (5), 1452–1463. Liu, E., Zhang, D., Liu, W., Liu, Y., 2009. Water consumption of the crops in the northwestern Shandong plain and comparisons of lysimeter and eddy covariance technique. Adv. Water Sci. 20, 190–196 (in Chinese with English abstract). Liu, C., Zhang, X., Zhang, Y., 2002. Determination of daily evaporation and evapotranspiration of winter wheat and maize by large-scale weighing lysimeter and micro-lysimeter. Agric. For. Meteorol. 111, 109–120. Makino, A., 2003. Rubisco and nitrogen relationships in rice: leaf photosynthesis and plant growth. Soil Sci. Plant Nutr. 49, 319–327. Maroco, J.P., Breia, E., Faria, T., 2002. Effects of long-term exposure to elevated CO2 and N fertilization on the development of photosynthetic capacity and biomass accumulation in Quercus suber L. Plant Cell Environ. 25, 105–113. Misson, L., Limousin, J.M., Rodrigues, R., Letts, M.G., 2010. Leaf physiological responses to extreme droughts in Mediterranean Quercus ilex forest. Plant Cell Environ. 33, 1898–1910. Mo, X.G., Liu, S.X., 2001. Simulating evapotranspiration and photosynthesis of winter wheat over the growing season. Agric. For. Meteorol. 109, 203–222. Mo, X.G., Liu, S.X., Lin, Z.H., Xu, Y., Xiang, Y., McVicar, T.R., 2005. Prediction of crop yield, water consumption and water use efficiency with a SVAT-crop growth

334

S. Hu et al. / Agricultural and Forest Meteorology 198–199 (2014) 320–334

model using remotely sensed data on the North China Plain. Ecol. Model. 183, 301–322. Mo, X.G., Liu, S.X., Lin, Z.H., 2012. Evaluation of an ecosystem model for a wheat–maize double cropping system over the North China Plain. Environ. Model. Softw. 32, 61–73. Moulin, S., Zurtia-Milla, R., Guerif, M., Baret, F., 2003. Characterizing the spatial and temporal variability of biophysical variables of a wheat crop using hyperspectral measurements. In: Proceedings of International Geoscience and Remote Sensing Symposium (IGARSS), pp. 2206–2208. Muller, J., Wernecke, P., Diepenbrock, W., 2005. LEAFC3-N: a nitrogen-sensitive extension of the CO2 and H2 O gas exchange model LEAFC3 parameterized and tested for winter wheat (Triticum aestivum L.). Ecol. Model. 183, 183–210. Norman, J.M., Becker, F., 1995. Terminology in thermal infrared remote sensing of natural surfaces. Agric. For. Meteorol. 77, 153–166. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2007. Numerical Recipes: The Art of Scientific Computing (3rd ed.), Section 10.2. Golden Section Search in One Dimension. Cambridge University Press, New York, ISBN 978-0-521-880688. Sage, R.F., 2002. Variation in the Kcat of Rubisco in C3 and C4 plants and some implications for photosynthetic performance at high and low temperature. J. Exp. Bot. 53 (369), 609–620. Schlerf, M., Atzberger, C., 2006. Inversion of a forest reflectance model to estimate structural canopy variables from hyperspectral remote sensing data. Remote Sens. Environ. 100 (3), 281–294. Sellers, P.J., Berry, J.A., Collatz, G.J., Field, C.B., Hall, F.G., 1992. Canopy reflectance, photosynthesis and transpiration. Part III: A re-analysis using improved leaf models and a new canopy integration scheme. Remote Sens. Environ. 42, 187–216. Sellers, P.J., Randall, D.A., Collatz, G.J., Berry, T.A., Field, C.B., 1996. A revised land surface parameterization (SIB2) for atmospheric GCMs. Part I: Model formulation. J. Clim. 9, 676–705. Shabanov, N.V., Knyazikhin, Y., Baret, F., Myneni, R.B., 2000. Stochastic modeling of radiation regime in discontinuous vegetation canopies. Remote Sens. Environ. 74 (1), 125–144.

Shuttleworth, W.J., Wallace, J.S., 1985. Evaporation from sparse crops – an energy combination theory. Quart. J. R. Meteorol. Soc. 111, 839–855. Tester, M., Davenport, R., 2003. Na+ tolerance and Na+ transport in higher plants. Ann. Bot. 91, 503–527. Tucker, C.J., 1979. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ. 8, 127–150. von Caemmerer, S., Farquhar, G., 1981. Some relationships between the biochemistry of photosynthesis and the gas exchange of leaves. Planta 153, 376– 387. Wan, Z.M., 2008. New refinements and validation of the MODIS land-surface temperature/emissivity products. Remote Sens. Environ. 112, 59–74. Wang, L., Koike, T., Yang, K., Yeh, P.J., 2009. Assessment of a distributed biosphere hydrological model against streamflow and MODIS land surface temperature in the upper Tone River Basin. J. Hydrol. 377, 21–34. Wang, Y.J., Tian, Y.H., Zhang, Y., E.I-Saleous, N., Knyazikhin, Y., Vermote, E., Myneni, R.B., 2001. Investigation of product accuracy as a function of input and model uncertainties – case study with SeaWiFS and MODIS LAI/FPAR algorithm. Remote Sens. Environ. 78 (3), 299–313. Wilson, K.B., Baldocchi, D.D., Hanson, P.J., 2001. Leaf age affects the seasonal pattern of photosynthetic capacity and net ecosystem exchange of carbon in a deciduous forest. Plant Cell Environ. 24, 571–583. Wullschleger, S.D., 1993. Biochemical limitation to carbon assimilation in C3 plants – a retrospective analysis of the A/CI curves from 109 species. J. Exp. Bot. 44 (262), 907–920. Xiong, W., 2009. The performance of CERES-wheat model in wheat planting areas and its uncertainties. J. Appl. Meteorol. Sci. 20 (1), 88–94 (in Chinese with English abstract). Xu, L., Baldocchi, D., 2003. Seasonal trends in photosynthetic parameters and stomatal conductance of blue oak (Quercus douglasii) under prolonged summer drought and high temperature. Tree Physiol. 23, 865–877. Zhao, C., Liu, L., Wang, J., Huang, W., Song, X., Li, C., Wang, Z., 2004. Methods and application of remote sensing to forecast wheat grain quality. In: Proceedings of International Geoscience and Remote Sensing Symposium (IGARSS), pp. 4008–4010.