Attributing regional trends of evapotranspiration and gross primary productivity with remote sensing: a case study in the North China Plain

Attributing changes in evapotranspiration (ET) and gross primary productivity (GPP) is crucial for impact and adaptation assessment of the agro-ecosystems to climate change. Simulations with the VIP model revealed that annual ET and GPP slightly increased from 1981 to 2013 over the North China Plain. The tendencies of both ET and GPP were upward in the spring season, while they were weak and downward in the summer season. A complete factor analysis illustrated that the relative contributions of climatic change, CO2 fertilization, and management to the ET (GPP) trend were 56 (−32) %, −28 (25) %, and 68 (108) %, respectively. The decline of global radiation resulted from deteriorated aerosol and air pollution was the principal cause of GPP decline in summer, while air warming intensified the water cycle and advanced the plant productivity in the spring season. Generally, agronomic improvements were the principal drivers of crop productivity enhancement.


Introduction
Terrestrial hydrological and carbon cycles are intimately coupled via transpiration and photosynthesis processes which are regulated by plant leaf stomata. Due to land use/cover changes, intensified agricultural management, and climatic change, terrestrial eco-hydrological processes have been noticeably shifted on multiple spatiotemporal scales (Tian et al., 2011;Douville et al., 2013); for example, prevailing irrigation and application of chemical fertilizers have raised soil moisture, evapotranspiration (ET), and crop productivity, etc. In some regions the effects of human activities are the same magnitude as, or even exceed, the impacts of global warming on the productions of agro-ecosystems (Haddeland et al., 2014). In the last decades, global consumptive water use and carbon fixation by terrestrial ecosystems have been demonstrated to slightly increase with more efficient water use, corresponding to changes in climatic factors and the fertilization effect of elevated atmospheric CO 2 concentration (Yan et al., 2013;Nayak et al., 2013). Consequently, spatiotemporal patterns of water and carbon fluxes on the regional scale are changing under global change (Zeng et al., 2014;Liu et al., 2012).
As ET is the major component of water budget in the water-limited basins, its long-term tendency has been taken as a critical indicator for diagnosing the intensification of the regional water cycle. The complementary relationship between actual and potential ET may reveal some clues to hydrological changes. Observations illustrated that potential evaporation rates (ET p ) (represented as pan evaporation) decreased in Europe, the US, China, India, and Australia in recent decades (Brutsaert, 2006;Katul et al., 2012), implying the decline of available energy and aerodynamics devoted to latent heat flux over the land surface. The climatic factors dominating ET p change are usually diverse. For example, over the North China Plain (NCP) the changes in ET p are mainly attributed to declines in global radiation and nearsurface wind speed (Tang et al., 2011;Song et al., 2009). However, in southern Turkey a noticeable decline in ET p was attributed to enhanced air humidity associated with the Published by Copernicus Publications on behalf of the European Geosciences Union. 296 X. Mo et al.: Attributing regional trends of evapotranspiration expansion of irrigation acreage and more water evaporated into the atmospheric boundary layer (Ozdogan and Salvucci, 2004). Burn and Hesch (2007) revealed that decreasing wind speed and raised water vapour deficit responded to the trend of ET p in the Canadian Prairies. On a large scale, precipitation is usually the principal factor determining actual ET change; for example, Qian et al. (2007) showed that an increase in ET in the Mississippi River basin followed precipitation propensity, while the effects of solar radiation and air temperature changes were minor.
Terrestrial eco-hydrological processes are driven by climate and modulated by human activities. Generally climate warming enhances atmospheric evaporative demand, while CO 2 fertilization stimulates photosynthesis and inhibits leaf stomatal conductance, leading to more biomass accumulation and higher water productivity (Field et al., 1995;Buckley and Mott, 2013). Simultaneously, land use change and land management also noticeably affect the ecosystem production and hydrological fluxes (Shi et al., 2011). Separating the contributions of climatic change, CO 2 enrichment, and human activities to the long-term trends of water and carbon cycles is critical for assessment of ecosystem responses and resilience to environmental changes. Some researchers have explored the relative contributions of climate change and vegetation dynamics to changes in global land surface evapotranspiration and river runoff (Betts et al., 2007;Piao et al., 2007;Alkama et al., 2010;Liu et al., 2012;Chen et al., 2014Chen et al., , 2015Banger et al., 2015), but the conclusions are still inconsistent. Climate change dominated the interannual variability of ET, while land use changes and agricultural practices and techniques exerted more discernable effects on the water cycle in the long term . However, Alkama et al. (2010) and Shi et al. (2011) demonstrate that climate change is the predominant driver of the global ET tendency in the 20th century. Commonly, the contributions of climate change to vegetation productivities on a large scale are quantified by the ecosystem models or statistical models. For example, Piao et al. (2015) documented that elevated atmospheric CO 2 and nitrogen deposition were the critical contributors to terrestrial greening over China; Baker et al. (2010) figured out that climate anomalies in springtime were the most frequent drivers of annual GPP variations in North America; Nayak et al. (2013) reported that climate change had a relatively small but significant control (15 %) on the trend of terrestrial net primary production (NPP) over India. In the crop ecosystems, contributions of climate change, cultivar renewal, and agronomic management to change in crop yield have been separated with the crop or statistical models (Yu et al., 2012;Song et al., 2014;Bai et al., 2016;Guo et al., 2014;. The impacts of climate change on crop yield may be positive or negative in different regions, depending on the tendencies of the dominant factors (Ewert et al., 2015).
As one of the granaries in China, the North China Plain (NCP) is confronted with challenges of agriculture sustain-ability due to global change and social and economic development. Thereby, it is crucial to understand the impacts of climate change on the productions of cropping systems. Over the plain, a winter wheat-summer maize double cropping system prevails, supported by irrigation, fertilizer, and agronomic techniques. In situ measurements, agricultural annals, and the regional remotely sensed vegetation index dataset all illustrate that both wheat and maize productivities have accrued remarkably since the 1970s (Mo et al., 2009;Yuan and Shen, 2013); correspondingly, the seasonal water consumption and water use efficiency are also slightly improved . The achievement of long-term increasing grain production is related to the active adoption of new varieties for stabilization, extending the length of the crop growth period, as well as agronomic technique advancement (Liu et al., 2010;Sacks and Kucharik, 2011). Currently, water amount for food production consists of 65 % of total water consumption here. Further, along with the gradually augmented domestic and industrial water requirement, groundwater in some areas of the plain has been overexploited, and the environmental water requirement is generally under deficit conditions (e.g. MWR, 2010). Faced with the rapid deteriorating agricultural environment, some critical issues are still unclear, such as what mechanisms drive the evolutions of eco-hydrological processes over the plain? What are the impacts of climate change on the cropping systems?
In this study, the VIP eco-hydrological dynamic model integrated with the NOAA-AVHRR remotely sensed normalized difference of vegetation index (NDVI) is employed to assess the spatiotemporal evolutions of ET and vegetation GPP over the NCP from 1981 to 2013. By numerical experiments with the VIP model and the factor separation method, the contributions of climate change, fertilization of atmospheric CO 2 enrichment, agronomic practices, and technological advancement to crop water consumption and productivity are then analysed and the relevant mechanisms are discussed.

Study region
The NCP is one of the country's granaries, extending from latitude 32 • 00 to 40 • 24 N and longitude 112 • 48 to 122 • 45 E (Fig. 1a, b). It is located in the eastern part of China, with an area of 40 × 10 4 km 2 , which is an alluvial plain developed by the intermittent flooding of the Huang, Huai, and Hai rivers, and 85 % is cultivated as farmland. The warm temperate climate varies gradually from subhumid in the southern to semi-arid in the northern parts. The annual precipitation ranges from 500 to 1000 mm, occurring irregularly among the seasons, and more than 70 % falls in summer. Soil moisture deficit happens widely during the spring and early summer period. Besides soy-bean/millet/sorghum/cotton, the double cropping system of winter wheat-summer maize prevails in the plain, where wheat and maize are the most common harvest crops in the summer and autumn seasons, respectively. Due to insufficient precipitation, the spring crops (such as wheat) usually need supplemental irrigation to form favourable production.

The VIP eco-hydrological model
The physical-process-based VIP (Vegetation Interface Processes) eco-hydrological model is designed to simulate the exchanges of energy, water, and carbon between the terrestrial ecosystem and the atmosphere ). In the model, ET is termed the summation of canopy transpiration and evaporation from the canopy intercept and the soil surface, computed separately with the Penman-Monteith equation. Transpiration and photosynthesis processes are coupled through the Ball-Berry relationship between leaf stomatal conductance and net assimilation rate. Regarding the carbon cycle aspect, leaf carbon fixations on sunlit and shaded leaves are predicted with the biochemical schemes for C3 (Farquhar et al., 1980) and C4 plants (Collatz et al., 1992). In the radiation budget scheme, shortwave radiation transfer in the canopy distinguishes the leaf spectral properties of visible and near-infrared radiation, as well as the fraction of direct beam and diffusive irradiance in global radiation. Precipitation throughfall-runoff generation over the land surface is calculated with a curve-number (CN) type equation on a daily scale, using the daily net precipitation and the moisture deficit of the upper soil layer. Simulation of soil water movement in the root zone is carried out with a discrete Richards equation in three layers. The crop and natural vegetation growth modules are also embedded in the model to simulate the biomass accumulation and carbon cycle.

Data
The VIP model input data include land use/cover, soil physical properties, and atmospheric forcing variables. The GIMMS AVHRR 15-day normalized difference of vegetation index time series (NDVI3g; values are scaled by 1000) (https://nex.nasa.gov/nex/projects/1349/wiki/general_ data_descriptionand_access/) (Pinzon and Tucker, 2014) from 1981 to 2013 is used to retrieve the vegetation leaf area index and other land surface characteristics. The land use classification originates from both Landsat TM images (www.resdc.cn) and MODIS remote sensing products, in which the farmland is classified as rice paddy and dryland. Soil textural data are available on a scale of 1 : 1 000 000 represented as fractions of sand, silt, and clay, by which the parameters of soil porosity (θ sat ) and saturated hydraulic conductivity (K ws , mm s −1 ) are estimated as in Bonan (1996). Daily climate variables (air temperature, water vapour pressure, wind speed, sunshine duration, and precipitation) recorded at 87 climatic stations (http://data.cma.cn) in and around the study area are available for generating the spatial atmospheric forces. The NDVI data are error-checked and the erroneous data are replaced by interpolation with the preceding and subsequent values according to the time series by the Savitzky-Golay (SG) filter (Savitzky and Golay, 1964), and then the daily values are derived with the Lagrange polynomial. Vegetation leaf area index (LAI) is retrieved from the NDVI with empirical relationships for different plant functional types.
The data used for model validation are field flux measurements with an eddy covariance technique at the Yucheng ( sites over the plain. The cropping systems at the Yucheng, Daxing, and Guantao sites are all rotation of winter wheatsummer maize, while land cover is dwarf shrub at the Miyun site. The eddy covariance data are processed with general procedures . GPP data are available only at the Yucheng site. In addition to eddy covariance fluxes, grain yield records of wheat and maize in county statistics are also used to verify the GPP predictions at the regional scale. The method of crop yield converted to GPP follows Xin et al. (2015).

Simulation setup
The model simulations were conducted at 8 km spatial resolution and a half-hour time step. The cropland practices are classified into wheat-maize and wheat-rice rotations. Atmospheric driving forces are interpolated from daily meteorological variables recorded at the climatic stations to grid cells with a gradient inverse distance square method (GIDS), which accounts for the effects of elevation, latitude, and longitude (Nalder and Wein, 1998). Estimated with sunshine duration in a linear relationship, the global radiation is subdivided into direct visible and near-infrared parts, as well as direct beam and diffusive components with Weiss and Norman (1985). The daily air temperature is extended to hourly values with a sinusoidal function based on the daily maximum and minimum temperatures (Campbell and Norman, 1998). During the winter wheat growing period, irrigation water is supplied when water storage in the root zone is below 60 % of the field capacity. Summer maize is set to be irrigated not more than one time in its growth period. The simulation is conducted with prescribed daily LAI series retrieved from remotely sensed NDVI series for eco-hydrological prediction from 1980 to 2013, in which the first year is taken as warming up.

Separation of climate change and management effects
By using a general function, f , the scalar fluxes (water vapour and carbon) between land surface and the atmosphere are determined by climate factors (M), fertilization of elevated atmospheric CO 2 (C), and agronomic management and technological advancement (in this study we assume the long-term trend of leaf area index (LAI) may represent the effects of human activities on crop and natural ecosystems. Human activities (also referred to as management practices) in the agro-ecosystem include renewals of cultivars, irrigation facility improvement, fertilizer application, soil quality amelioration, etc.), namely, The changes in f contributed by a single factor (expressed as f i ) and its interaction with another factor (expressed as f ij ) can be decomposed by the Taylor expansion as where x i and x j (i = j ) represent M, C, and LAI, respectively. The factor separation methodology from Stein and Alpert (1993) and Alkama et al. (2010) is used to categorize the contributions of climate change, CO 2 fertilization, and LAI, as well as their interactions to long-term trends of ET and GPP. Similar to Alkama et al. (2010), the total effect, f 123 , is expressed as with where f 1 , f 2 , and f 3 are the direct contributions of climate change, atmospheric CO 2 enrichment fertilization, and agronomic management, respectively; f 12 is the contribution of interactions of climate change and CO 2 enrichment; f 13 is the contribution of interactions between climate change and agronomic management; f 23 is the contribution of interactions between CO 2 fertilization and agronomic management; f 123 is the contribution of interactions between climate change, CO 2 fertilization, and agronomic management. Seven numerical experiments designed to fully distinguish the contributions of climate change, CO 2 fertilization, and agronomic management are conducted by the VIP model over the NCP from 1980 to 2013. The experiments are as follows: 1. f 123 ("all factors"): current climate, CO 2 , and LAI spatiotemporal pattern;  5. f 12 ("climate change and CO 2 fertilization effects"): LAI pattern is fixed, but the current climate and CO 2 concentration are used; 6. f 13 ("climate change and management effects"): CO 2 concentration is fixed, but the current climate and LAI pattern are used; and 7. f 23 ("CO 2 fertilization and management effects"): climate is fixed at 1981, but current CO 2 and LAI are used.
The trends of annual ET and GPP in the above experiments are calculated, and then, according to Eqs. (4) to (7), the contributions of climate change and CO 2 fertilization and management practices to ET and GPP long-term trends are then separated.

.1 Validation with eddy covariance measurements
The VIP model is used to simulate the hydrological, energy partitioning, and crop growth processes at the four sites of eddy flux measurement. Here, eddy covariance measurements of daily ET and GPP are employed to verify the model predictions (ET is available at all the sites, but GPP is only available at one site). The land surface characteristics surrounding the flux mast are relatively homogeneous, ensuring the footprint for flux measurement. The meteorological information measured at each site is used to drive the VIP model. It is shown that the agreements are quite satisfactory for both ET and GPP (Figs. 2 and 3). In total, there are 9-year daily ET data and 3-year daily GPP data for comparison with the model simulations. The coefficients of determination (R 2 ) are above 0.76 for all the sites. On an annual scale, the relative absolute biases of predicted ET ranged from 1.5 to 12.6 % in the 9-year dataset, and biases of GPP are from 2.0 to 8.8 % in 3-year data. Therefore, the model performance is quite good and reliable for predictions of vegetation/crop productivity and water consumption. The biases may have stemmed from both measurements and model parameter uncertainty. Mo et al. (2012) showed that canopy leaf area index (LAI) and photosynthetic capacity (carboxylation rate for C3 crops and photon quantum use efficiency for C4 crops) were the parameters most sensitive to the model efficiency. Here, taking the Yucheng site as an example, annual ET and GPP may increase by 1.6 % (2.6 %) and 3.0 % (15.9 %), respectively, as LAI (photosynthesis capacity) is increased by 20 %.

Validation with the statistical yield records
The simulated GPP is also validated with the staple crop grain yield statistics at county level. The crop grain yield per hectare is converted to the equivalent GPP per square metre. As shown in Fig. 4 (years 2000 and 2005 are used), the Figure 5. Spatial average trends of the growing season NDVI at annual (a) and monthly (b) scales (significance levels: * is p < 0.05; * * is p < 0.01; the NDVI is scaled by 1000).
agreement is satisfactory, with coefficients of determination (R 2 ) of 0.43 and 0.51 (p < 0.001), respectively. It is noted that the model overestimated low-yield but underestimated high-yield zones. The errors may have resulted from two aspects, namely, biases in the LAI-NDVI empirical relationship and the harvest index in low-and high-yield fields. There are remarkable spatial variations in crop yields that resulted from diverse climates, soil physical texture, chemical composition, and management practices. In the simulation, it is found that the spatiotemporal evolution of greenness is a crucial factor regulating the yield patterns. Greenness represented by the vegetation index is an appropriate indicator of crop productivity under environmental stresses (Hu et al., 2014). In the areas with a high vegetation index and favourable irrigation facilities, the yield losses may be caused by heat waves or pest infections in the mature stage.

Trends of climate, crop productivity, and ET
Changes in climate variables and agro-ecosystem management are the dominant driving forces of the evolution of regional eco-hydrological processes. Intra-seasonal variations of climatic variables may exert different impacts on the crop water consumption and carbon assimilation. In the last 3 decades, air temperature has been rising, but sunshine duration and wind speed are decreasing significantly over the plain, associated with global climate change, aerosol, and air pollution exacerbation. Soil amelioration, genetic improvement, irrigation facility constructions, and application of chemical synthesis fertilizer are considered to be the principal drivers that have propelled the crop productivity close to the attainable level (Yu et al., 2012;Lobell and Burke, 2010).

Changes in climatic variables
Grid averages of the climatic variables were analysed over the North China Plain from 1980 to 2013. Nevertheless, inhomogeneous distributions of the climatic variables and the spatially averaged trends were clear (Table 1). On an an-nual scale, global radiation, air temperature (especially minimum temperature), and wind speed changed significantly (p < 0.01). On a monthly scale, radiation declined in all the months except March, but only trends in June to September were significant (p < 0.01); a significant increase in air temperature occurred in spring (February and March) and early summer (May to July); wind speed decreased significantly (p < 0.01) in all the months except August. However, no significant trends were detected for both precipitation and water vapour pressure throughout each month. As a consequence, water vapour pressure deficit was exaggerated along with the rise in air temperature, which was expected to intensify the atmospheric water vapour demand and offset the negative effects of declining radiation and wind speed on potential evaporation. These changes in climatic variables have exerted remarkable impacts on the crop phenological stages, water consumption, and productivity during the last 3 decades over the North China Plain (Liu et al., 2010).

Changes in greenness and GPP
Here the remotely sensed NDVI is expressed as vegetation greenness. Averaged over the growth period (from March to October), vegetation greenness significantly increased from the 1980s, with a trend of 0.64 yr −1 (p < 0.001) and a coefficient of variation (CV) of 2.4 % (Fig. 5a, b). The maximum inter-annual variation of greenness was 5.4 % between 2 consecutive dry and wet years (1989 and 1990), with a 280 mm difference in annual precipitation. The distinct differences in greenness occurred in June to December. On an annual scale, greenness was weakly related to precipitation; however, in the growing season, greenness was noticeably correlated with monthly precipitation in April, May, June, September, and November (r = 0.29-0.53, p < 0.1). The reason is that the monthly rainfall is generally lower than the atmospheric evaporative demand in spring, and the water deficit to transpiration generally stresses the plant growth. Unlike precipitation, greenness anomalies are positively correlated with the detrended air temperature (r 2 = 0.16, p < 0.05), implying that recent climate warming has stimulated vegetation growth by extending the growing stage and by pushing photosynthesis in water non-limited regions (Mao et al., 2012;Nemani et al., 2003).
Regional greenness trends were remarkably diverse (Fig. 6a, b). Except for climate change, human activities also exerted critical impacts on the terrestrial greenness variations. In the low plain of Hebei province, saline-alkali land amelioration and irrigation facility improvement contributed greatly to the vegetation density and greenness enhancement in the 1980s to the 1990s. In addition, atmospheric nitrogen deposition is also regarded as a positive driver of the vegetation improvement, since the nitrogen deposition increased on average by 25 % from the 1990s to the 2000s in North China (Jia et al., 2014;Piao et al., 2015). Spatially, greenness increased over 91.3 % of the NCP, in which the most distinc- Table 1. Inter-annual trends of monthly sunshine duration (Sun), precipitation (P ), air temperature (T ), relative humidity (RH), and wind speed (U ). tive grids were distributed in the southern parts and the belt along the Yellow River channel, where water supply was usually sufficient. In the northern part of the plain, however, the tendencies of greenness in a number of grids decreased significantly at the p = 0.1 level, which were assumed to have resulted from less irrigation supply to crops in springtime and rapid expansion of built-up occupations around cities and towns, such as the Beijing and Tianjin metropoles. The spatially averaged GPP was 1913 ± 584 gC m −2 yr −1 , with a CV of 6.8 % predicted by the VIP model from 1981 to 2013 showing noticeable spatial variability (Fig. 7). Low crop productivity resulted from fields with saline-alkali soil in the low lands near the coast of the Bohai Sea, where almost no favourable water was available for irrigation purposes in springtime. On average, the increasing trend of GPP was significant, with a slope of 8.2 gC m −2 yr −2 (r = 0.60, p < 0.01). It was noticed that the average annual GPP increased steadily from the 1980s to the 2000s, compounded by decadal variations of climate and elevated atmospheric CO 2 , as well as the improvement of agricultural practices and techniques. It was revealed that trends of annual GPP were positive over 87.9 % of the study region. As shown in Fig. 7, the obvious increasing trends were located in the mid and southern areas, while most of the decreasing trends occurred in the eastern and northern parts, where water for irrigation was considerably reduced in the spring season because of the competing demands of domestic and industrial water use.
On a monthly scale, GPP increased in all the months except for July, August, and September (Fig. 8). The positive trends were contributed principally by the summer harvest crops (wheat as the main crop), while the negative trends were mainly contributed by the autumn harvest crops (maize as the major crop). Regressive analysis showed that the downward trends of GPP in the summer season resulted from the significant declines in monthly sunshine duration and shortwave radiation (r = 0.38 to 0.57 from June to August, p < 0.05).

Changes in ET
Water loss from the vegetation surface as ET is directly regulated by atmospheric vapour demand and leaf stoma physiological functioning (Buckley and Mott, 2013). Inter-annual variation of ET is controlled by climate variability/change and agronomic managements. Generally potential ET (ET p ) is used to represent the available energy for water vaporization on land surface. As shown in Fig. 9a, ET p was slightly decreasing, which resulted from offsetting between the negative effects of reduced global radiation and declining wind speed and the positive effect of increasing water vapour deficit (Song et al., 2009); simultaneously, actual ET was predicted to be slightly increasing (p < 0.05), consistent with the enhancement of greenness. It was noticed that the evolutions of potential and actual ET coincided with the hypothesis of complementary relationship. On a monthly scale, ET significantly increased from February to April, but it decreased in August (Fig. 9b). This implies that climate warming may be beneficial to spring crops by advanced recovering date from dormancy, whereas a decline of net radiation (R n ) (especially in August, significance level p < 0.001) may lead to the downward tendency of ET rates in summer.
Over the whole plain, spatially averaged actual ET and transpiration were 627 ± 162 mm yr −1 (about 92 % of annual precipitation) and 416 ± 129 mm yr −1 (about 67 % of ET), respectively. The trend of annual ET (p < 0.1) with CV of 0.05 was 0.88 mm yr −2 from 1981 to 2013, which was less significant than that of the NDVI (p < 0.01). Decadal ET amounts in the 1980s, 1990s, and 2000s were 610, 626, and 640 mm, respectively, corresponding to the slightly rising trend of precipitation. It is found that GPP increased with a higher significance level than that of ET, implying the enhancement of water productivity in the plain. Spatially, the trends of ET were positive over 86.0 % of the study region, mainly occurring in the mid and southern parts, while negative trends were mainly found in the northern part (Fig. 10a,  b), which was consistent with the pattern of GPP tendencies. By using a water balance model, Zeng et al. (2014) also presented an increasing trend of ET over the North China Plain from 1982 to 2009.

Contributions of climate change, atmospheric CO 2
fertilization, and agronomic management to changes in ET and GPP

Spatial patterns of the contributions
The contributions from climate change, atmosphere CO 2 enrichment fertilization, and agronomic management illus-  trated considerable spatial heterogeneity for both ET and GPP (Fig. 11). Over the whole plain, climate change exerted a positive impact on water vapour exchange from the land surface to the atmosphere (f 1 ), especially in the eastern hilly part where precipitation increased slightly. As is generally known, air CO 2 enrichment stimulates the crop leaf stomatal closing and then reduces transpiration, but its fertilization effect enhances the photosynthetic rate and water use efficiency (Buckley and Mott, 2013). Descriptions of the separate effects are presented as follows.
The climate change has intensified the hydrological cycle, resulting in a 0 to 4 mm increment of ET per year. The effect of climate change was much stronger in the mid to east-ern zones with high crop productivities, contributed mainly by air temperature increase. The contribution of CO 2 enrichment to ET is negative in most areas, ranging from 0 to −1 mm per year. The attributions of agronomic practices and technological advancement represented by LAI increase are somewhat complex, namely a remarkable increase ranging from 0 to 6 mm per year in the mid-western areas where irrigation facilities and soil conditions have been ameliorated greatly in the recent decades through land consolidation and de-salinization. Renewal of cultivars and improved agronomic practices also contributed to ET intensifying . In contrast, ET in the grids with expansions of built-ups appeared to have negative trends. Annually, the contribution of climate change to GPP is negative, ranging from 0 to −12 gC m −2 per year. Low rates occurred in the southwestern and northern parts. Air warming and heat waves, declines in precipitation, and global radiation are the main causes of crop production reduction (Lobell et al., 2011;Guo et al., 2014). In addition to the spatial variability, climate change effects are associated with the relevant land use/cover and cropping systems. In the hilly (western and mid) and eastern coast areas, negative effects were slight, where air warming and air pollution were relatively weak. CO 2 enrichment effects were positive over the whole plain, ranging from 0 to 6 gC m −2 per year. It was noticed that the higher effects were associated with higher cropland with favourable irrigation and high productivity, while the lower rate was related to low-productivity croplands and natural vegetation communities. Similarly, the effects of human activities on ET were positive in the mid to western areas, ranging from 0 to 35 gC m −2 per year, associated with croplands of high productivity. The negative effects mainly occurred in the eastern and northern parts, where there was remarkable expansion of urban and dwelling buildings in the study period.

Regionally averaged contributions
Regarding the aspect of regional average, some characteristics of the contributions to water and carbon assimilation are revealed. As shown in Fig. 12a, the contributions of climatic variable change (f 1 ), elevated atmospheric CO 2 concentration (f 2 ), and agronomic management (represented by leaf area index (LAI) increment) (f 3 ), and their interactions with the long-term trend of ET were positive, while the contribution of elevated atmospheric CO 2 was negative in the last 3 decades. It was shown that the contribution of climate change was less than that of agronomic improvement. The relative direct contributions of climatic change, CO 2 fertilization, and agronomic management and technological advancement to the long-term trend of evapotranspiration were 56, −28, and 68 %, respectively. Compared with the contributions of direct effects, the relative contributions by their in-teractions were low (the cumulative effect of f 12 , f 13 , f 23 , and f 123 was only about 4 %). Although the global radiation reaching the ground was diminished by a higher aerosol concentration and deteriorated pollution in the atmosphere (Che et al., 2005), its negative effect on terrestrial ET was offset by the positive effects of air warming and its related higher vapour pressure deficit (VPD) on ET on an annual scale. Reduction of transpiration by enriched atmospheric CO 2 caused by closure of plant leaf stomata at high CO 2 concentrations for both C3 and C4 plants may mediate the extra water demand by air warming. The dominant contribution to ET changes was from the renewal of cultivars and improvement of agricultural techniques and management. In the study period, agronomic management has greatly improved, including the establishment of irrigation facility, prevalent uses of chemical fertilizers, and pesticides. For example, irrigated area in the northern part (mainly the Hebei Plain) has increased by 2.5 times, and chemical synthetic fertilizer input has increased about 4 times; as a consequence, crop grain production enhanced about 2 times from the 1980s to the 2000s . Climate change and agronomic management improvement (irrigation practice, synthetic fertilizer supply, and new cultivar adoption) are the main contributors to ET intensifying over the plain.
As shown in Fig. 12b, the enriched atmospheric CO 2 fertilization and agronomic management improvement presented a positive contribution to the GPP trend during the study period. It was somewhat against expectations that the contribution of climate change to GPP was negative on an annual scale. The relative contributions of climate change, CO 2 fertilization, and management to the vegetation GPP enhancement were −32, 25, and 103 %, respectively. It was confirmed that the improvement of agricultural management was the dominant driver of GPP increase. The positive effects on GPP were associated with human activities and natural factors, such as input of synthetic fertilizers and atmospheric nitrogen deposition, irrigation, and other agronomic technology improvement, as well as fertilization of enriched atmospheric CO 2 . The negative contribution by climate change mainly happened in summertime (Fig. 8). Since there was less benefit of CO 2 enrichment to summer maize (C4 type), the reduced maize productivity due to global radiation decline was not fully offset. Some studies, such as Piao et al. (2015), also reported that climate change exerted a negative impact on the vegetation greening phenology, and Liu et al. (2010) attributed the reduction of crop productivity to shortening of vegetative growth length under climate warming. As shown in Fig. 4b, it was illustrated by NDVI time series that greenness in the summer season was quite stable; however, it significantly increased in spring and autumn, indicating that climate warming was beneficial for crop growing in the cool seasons. Additionally, carbon assimilated by summer crops was larger than that by spring crops. Thereby, the sign of the annual GPP trend was determined by its trend in the summer season. As current air temperature increase Figure 9. Inter-annual trends of potential and actual annual ET (a) and trends of monthly ET and net radiation (R n ) (b) from 1981 to 2013 (significance levels: * is p < 0.05; * * is p < 0.01). was not yet so detrimental to maize growth, the decline of downward shortwave radiation was considered to be responsible for GPP decline.

Effects of climatic variables on monthly ET and GPP trends
To attribute the responses of cropping systems to the trends of single climatic variables, the VIP model is used to diagnose the effects of climate change on ET and GPP at the Beijing meteorological observation site. The contributions of a single variable to the trends of ET and GPP are expressed by their differences simulated with the current and de-trended variables, respectively. Here only the climatic variables of radiation, air temperature, and wind speed were linearly detrended on a monthly scale, since no significant trends of precipitation and specific humidity are detected. As shown in Table 2, while the global radiation was de-trended, the negative correlation coefficient of monthly GPP with time was reversed from negative to positive in springtime, and from significant (r < −0.6, p < 0.01) to insignificant levels (r < −0.15, p > 0.01) in July and August. It was affirmed that the decline in global radiation was the dominant factor for reduction of crop GPP in the summer period (June to August), but in the autumn season the changes in radiation, temperature, and wind speed were all responsible for GPP changes. From Table 2, it could be deduced that the effect of temperature rising on crop productivity was positive and significant in spring (March and April) and autumn (September and October), whereas its effect was relatively weak in summer (May to August). It was noticed that the effect of radiation change was quite weak in March, when no significant trend of shortwave radiation was detected. In the spring season, sunshine duration increased from the 1980s to the 2010s (Wang and Yang, 2014). In June, except for global radiation, changes in temperature and precipitation have contributed to GPP increasing. Additionally, the fertilizing effect Figure 11. Contributions of climate changes, CO 2 enrichment, and management to ET (a, b, c) and GPP (d, e, f), respectively (a and d are for climate, b and e for CO 2 , and c and f for management). of enriched atmospheric CO 2 on C3 crops is a non-ignorable factor in GPP increasing. However, maize as the C4 crop does not benefit much from atmospheric CO 2 enrichment. It is suggested that new cultivars with higher light use efficiency should be adopted to sustain the maize productivity under declined global radiation condition resulting from exacerbating aerosol concentration and air pollution. Comparatively, the effect of climatic change on ET was less significant than that of vegetation GPP. The model sim-ulations showed that ET enhanced by air temperature rising mainly occurred in August to October, while the negative effect of solar radiation decreasing was detected from June to September in the maize growing period.

Discussion
Our simulations suggested that annual ET and vegetation GPP increased over the North China Plain from 1981 to Table 2. Monthly Pearson correlation coefficients of GPP trends (r _ ALL: all variables are not de-trended; r_R: radiation is de-trended; r_T : air temperature is de-trended; r_R_T _U : radiation (R), temperature (T ), and wind speed (U ) are all de-trended).

Jan
Feb 2013. Climate change contributed positively to ET intensification, but negatively to GPP enhancement. Agronomic management and technological advancement are the dominant factors in promoting GPP increase. The use of remote sensing NDVI series has greatly improved the reliability of vegetation water consumption and productivity prediction at spatial and temporal scales, even if there were uncertainties in vegetation characteristics retrieved from the NDVI dataset.
The results in this study were supported and consistent with the most relevant studies on field and regional scales.

Is the trend of ET upward or downward over the NCP?
Although the crop productivities are steadily increasing, whether the actual ET over the NCP is increasing or decreasing during the last 3 decades is controversial from the reported literatures. By using the complementary relationship models (Brutsaert and Stricker, 1979), actual ET and potential ET both decreased (Cao et al., 2014;Gao et al., 2012). However, ET increased, predicted by the process-based VIP model from 1981 to 2013, which was consistent with the increasing trend of terrestrial greenness (S. . Yuan and Shen (2013) found that in the northern part of the NCP (Hebei Province), ET was positively correlated with recorded crop grain yield, and agricultural water use has increased. Field measurements under well-watered fields also showed that seasonal ET rates of both winter wheat and summer maize increased . Brutsaert (2006) acknowledged that decreasing pan evaporation was evidence of increasing terrestrial evaporation. As is generally known, the sign of ET change should be the same as that of vegetation greenness. Over the NCP a positive trend of ET was more believable, in the light of a significantly increasing NDVI over the growing season, especially in spring. The positive effects of warming with higher water vapour deficit on ET might be offset by the negative effect of declining solar radiation and wind speed on potential evaporation. In view of the complementary relationship hypothesis (Hobbins et al., 2001), alteration of available energy partitioned into latent heat flux (or ET) is dominated by the atmospheric water vapour deficit, namely, while more vapour is evaporated into the atmosphere boundary layer, its water vapour deficit is then accordingly relaxed, resulting in a lower rate of ET p . However, while declining global radiation is the domi-nant factor in the ET trend, actual ET (ET a ), wet surface ET (ET w ), and potential ET (ET p ) trace the trend of available energy (net radiation). In this study case, net radiation declined at a rate of −5.58 MJ yr −1 (r = 0.56, p < 0.01) over the NCP.
As the trends of both ET p and ET w were dominated by the radiation trend, the ET a estimated from the complementary relationship definitely followed the negative trend of radiation, because the positive trend of aerodynamic evaporation was weak as a tradeoff of the positive effect of rising water vapour deficit and the negative effect of decreasing wind speed. However, the declining ET trend which resulted from the reduced radiation has actually been reversed by increasing green leaf area, which would reduce land surface albedo and temperature, etc. Consequently, ET and GPP showed a slight increase. This study case also confirmed limitations of the complementary relationship for diagnosing the evaporation trend under the condition of declining radiation.

Effects of climate change and CO 2 on the cropping systems
Without adaptation measures, climate change is shown to exert negative effects on the productivity of cropping systems in the NCP during the last 3 decades (Mo et al., 2013;Liu et al, 2010). Changes in individual climate variables affected specific cropping systems differently, which was associated with crop type and growing season. Climate warming in the winter and spring seasons was beneficial to vegetative growth of winter wheat (Mo et al., 2013). Although air warming has shortened the growing length, the autonomously adopted cultivars with higher thermal requirements usually maintained the crop growth length and accumulated more photosynthesis products, which may outweigh the extra respiration consumption under warmer climates . So far, the effects of global warming on wheat production have been inhomogeneous in the North China Plain, which were positive in the northern part but negative in the southern part . The reasons are that during the wheat growth period air temperature is still below the favourable conditions in the high-latitude part of the plain; therefore, recent global warming is benign to the wheat growth, whereas air warming may be detrimental in the southern part, especially for rainfed wheat (Xiao and Tao, 2014). However, dominated by summer monsoon in the North China Plain, the climate is hot in the summer maize growth period. Due to maize being a tropically originated species with a high thermal requirement, it can tolerate relatively high air temperature. Our study showed that maize has not suffered noticeably from air warming in the recent decades, confirmed also by Xiong et al. (2012). However, Guo et al. (2014) reported that effect of air warming on maize was adverse with an Agro-Ecological Zones model, and the decreased daily temperature range (DTR) may be detrimental to crop yields. Nevertheless, adaptation measures may boost the production, such as harvest time delay  or planting date advancement (Sacks and Kucharick, 2011).
As atmospheric aerosols from industrial production and combustion have increased, global radiation has declined in many parts of the world, in which direct components decreased but diffuse components increased, so-called "global dimming" (Liepert, 2002;Ren et al., 2013). The decline of global radiation has resulted in less pan evaporation and carbon assimilation in crop and natural vegetation communities (Xiong et al., 2012;Xiao and Tao, 2014); nevertheless, plant canopies can use the diffuse radiation with higher efficiency than direct beams (Gu et al., 2002). In springtime, while the atmospheric circulation is shifting from continental to ocean monsoon in East Asia, the wind speed is relative high, and consequently air pollution and aerosols are usually low; therefore, global radiation is not reduced obviously in the wheat growth period (Table 1), illustrating that air warming and precipitation variability other than radiation decrease are the principal climate factors contributing to the tendency of wheat production. In contrast, global radiation declines significantly in the summer season in the North China Plain (Table 1); as a result, productivities of autumn harvest crops such as maize are mainly affected. For example, Guo et al. (2014) reported that maize potential productivity was reduced by 20 kg hm −2 due to global radiation decline in the last decades over China.
During the study period, atmospheric CO 2 concentration increased from 340 to 396 ppm, which contributed to enhancement of crop productivity. However, most studies with statistical analysis models neglected the contribution of CO 2 fertilization (e.g. Lobell and Burke, 2010;Song et al., 2014). As confirmed by FACE experiments, elevated atmospheric CO 2 concentration is accelerating plant photosynthesis and reducing transpiration, whose fertilizer effect is 0.065 % per ppm increment for C3 plants (Field et al., 1995;Long et al., 2006;Ainsworth et al., 2008). In this study, the contribution of CO 2 to GPP is positive, but negative for ET. The positive CO 2 effect on GPP almost compensates for the negative effects of climatic variable changes. However, it should be borne in mind, as Ainsworth et al. (2008) pointed out, that the CO 2 fertilizer effect may be overestimated by the process-based crop/ecosystem models.

Effects of agronomic practice, technique advancement and other factors
The remotely sensed NDVI is an excellent indicator of longterm changes of vegetation covers. Here we assumed that climate change did not modify the tendencies of vegetation covers, but dominated the inter-annual variation of vegetation density. Renewal of crop cultivars, applications of synthetic fertilizer and irrigation, as well as conservancy tillage and nitrogen deposition all contribute to the crop and/or natural productivity improvements (Yu et al., 2012;Bai et al., 2016;Piao et al., 2015). National statistical records of grain yields on county scale showed a rapid increase from 1980 to 1990 and a moderate increase in the 2000s. Enhancement of crop yields mainly stemmed from more biomass accumulation and a higher harvest index than previous varieties . In our simulations the upward trend of GPP was more significant than that of ET, which was consistent with the increasing trend of the cumulative NDVI. Agronomic practices and technology advancement contributed 103 % of GPP changes in our study. By using crop models, Yu et al. (2012) and Song et al. (2014) reported, respectively, the relative contributions of 92 and 62 % by agronomic management and renewal of cultivars for rice, and Guo et al. (2014) showed that 99.6 to 141.6 % of maize yield increases was contributed by technological advancement in China since the 1980s. The previous studies showed that, if no adoption measures were taken, climate change generally contributed negatively to crop productivities in the mid-latitude areas, but the negative effects were usually compensated for by genetic improvements, applications of fertilizer and irrigation, pest and weed control, as well as CO 2 and nitrogen deposition effects (Liu et al., 2010;Lobell et al., 2011;Guo et al., 2014;Bai et al., 2015). Under warming climate conditions, it is expected that water requirement by crops and natural plants will increase, but the intensified ET may be limited by insufficient soil moisture availability. Therefore, the sustainability of crop productions greatly depends on the improvement of agronomic management and technological advancement on variety breeding.
Climate change and human activities have greatly altered the hydrological regime and crop productivity in the North China Plain, with warm temperate climate during the recent 3 decades. The VIP ecological model integrated with the NOAA-AVHRR NDVI data series predicted that spatial average annual actual ET weakly increased, while vegetation primary productivity (GPP) significantly increased (p < 0.01) from 1981 to 2013, consistent with the remotely sensed NDVI trend. The increases in actual ET and GPP mainly occurred in spring, while ET and GPP were obviously decreasing in August owing to global radiation diminishing. Climate change, elevated atmospheric CO 2 fertilization, and agronomic management (new cultivars, irrigation, and chemical fertilizer) all contribute to the inter-annual trends of ET and crop GPP. The relative direct contributions of climatic change, CO 2 fertilization, and agronomic management to ET increase are 56, −28, and 68 %, while the contributions to GPP tendency are −32, 25, and 103 %, respectively. Air warming intensifies the crop water requirement and enhances the production of crops harvested in summer. The decline of global radiation resulting from exaggerated aerosol concentration and air pollution is considered to be the main cause of GPP reduction in August. The study confirms the necessary for imminent control of air pollution and aerosol to sustain the productivity of the agricultural system.

Data availability
The VIP model simulation and analysed data products are available from the corresponding author. However, the measurement data of eddy covariance and yield census data are from some specific data providers with agreements; these data have limited accessibility.