GLOBALLY INCREASED CROP GROWTH AND CROPPING INTENSITY FROM THE LONG-TERM SATELLITE-BASED OBSERVATIONS

Understanding the spatiotemporal change trend of global crop growth and multiple cropping system under climate change scenarios is a critical requirement for supporting the food security issue that maintains the function of human society. Many studies have predicted the effects of climate changes on crop production using a combination of filed studies and models, but there has been limited evidence relating decadal-scale climate change to global crop growth and the spatiotemporal distribution of multiple cropping system. Using long-term satellite-derived Normalized Difference Vegetation Index (NDVI) and observed climate data from 1982 to 2012, we investigated the crop growth trend, spatiotemporal pattern trend of agricultural cropping intensity, and their potential correlations with respect to the climate change drivers at a global scale. Results show that 82.97% of global cropland maximum NDVI witnesses an increased trend while 17.03% of that shows a decreased trend over the past three decades. The spatial distribution of multiple cropping system is observed to expand from lower latitude to higher latitude, and the increased cropping intensity is also witnessed globally. In terms of regional major crop zones, results show that all nine selected zones have an obvious upward trend of crop maximum NDVI (p < 0.001), and as for climatic drivers, the gradual temperature and precipitation changes have had a measurable impact on the crop growth trend.


NTRODUCTION
Propelled by a 2.3-billion global population growth and higher per capita incomes anticipated through the mid-21 st century, global demand for agricultural crops is increasing and may continue to witness an upward trend for decades (Godfray et al. 2010).Food provision serves a prerequisite for the function of human society, and cropland where food and feed are grown is the central, limiting resource for food production (Kastner et al. 2012).Thus, the crop growth situation with regard to the global cropland areas is sensitively correlated to the food security issue including the four key dimensions of food supplies: availability, stability, access, and utilization (Food and Agriculture Organization (FAO), www.fao.org).Besides, crop production gains can also come from intensification, such as increased cropping intensity, and shorter fallow periods.The spatiotemporal pattern change of agricultural cropping intensity has a direct impact on the crop production and its related food security issue.All these emergent and potential issues point to the need for comprehensive monitoring and inventorying efforts on the global spatiotemporal change trend of long-term crop growth and different cropping system distribution.
Global crop monitoring is a long-term, large-scale complicated scientific project.Fortunately, satellite remote sensing has greatly facilitated the mapping and monitoring of croplands by providing spatially explicit and temporally continuous observations.Over the past decades, a number of studies have utilized remotely sensed data to extract cropland extents, quantify crop types, estimate crop yields, and monitor crop growth trend (Benedetti and Rossini 1993;Biradar et al. 2009;Dong et al. 2015;Fritz et al. 2010;Pittman et al. 2010;Seelan et al. 2003;Siebert et al. 2010).However, the majority of previous studies focuses on either regional/continental scales, or seasonal/short-term period.The long-term crop growth trend at a global scale is still limited, and the spatiotemporal pattern of global agriculture intensification is poorly characterized and understood.Being one of the top threats for the Earth in the 21 st century, the magnitude, rate, and pattern of climate change also greatly impacts on agricultural productivity.Crop growth is affected by biophysically by meteorological variables, including rising temperatures, changing precipitation regimes, and increased atmospheric carbon dioxide levels (Parry et al. 2004).Biophysical effects of climate change on agricultural production will be positive in some agricultural systems and regions, and negative in others, and these effects will vary through spatial and temporal difference (Parry et al. 2004).
To address the aforementioned issues, here we employed longterm Global Inventory Modelling and Mapping Studies (GIMMS) dataset from 1982 to 2012 to provide a global crop monitoring with the following three major objectives in this article: (i) examine global crop growth trend over the past three decades; (ii) investigate spatiotemporal pattern change of multiple cropping system; and (iii) attribute the major drivers of crop growth trend within climate change scenarios.

Datasets
The Normalized Difference Vegetation Index (NDVI), defined as the ratio of the difference between near-infrared and red visible reflectance to their sum, is a remotely sensed vegetation index widely used to measure vegetation greenness (Myneni et al. 1997;Tucker 1979).Here we used the GIMMS third generation biweekly NDVI dataset derived from AVHRR sensors (NDVI3g) with a spatial resolution of 8 km from 1982 to 2012 (Tucker et al. 2005).This GIMMS NDVI dataset has been corrected to reduce the deleterious variation arising from the calibration, viewing geometry and volcanic eruptions (Piao et al. 2011;Tucker et al. 2005), and has been widely used in terrestrial net primary production (NPP) estimation, vegetation phenology monitoring, and continental scale trends associated with environmental changes.The MODIS monthly global vegetation product dataset (MOD13C2) (Solano 2010) at a spatial resolution of ~5.6 km, spanning from 2000 to 2012 was also integrated.In order to be spatially and temporally comparable to the GIMMS and MODIS dataset, we acquired the AVHRR GLCF map at a spatial resolution of ~8 km (De Fries et al. 1998), which is utilized to mosaic cropland type from the time-series GIMMS dataset.In this way, we can greatly eliminate the impacts from other land cover types such as forests and meadows.
The CRU TS 3.0 climate dataset including monthly temperature and precipitation dataset that spanning from 1982 to 2012 was obtained from the Climate Research Unit (CRU) at the University of East Anglia (Mitchell and Jones 2005).This gridded dataset with a spatial resolution 0.5° x 0.5°, was based on climate observations from more than 4000 meteorological stations (Mitchell and Jones 2005;Wang et al. 2011).
The Palmer Drought Severity Index (PDSI), one of the most commonly used drought indices, was adopted to indicate spatiotemporal variations of drought.The monthly PDSI dataset produced by Dai et al. (Dai et al. 2004) with a spatial resolution of 2.5° x 2.5° was used in this study, which also covers the entire study period .
Nine major crop zones were defined for regional analysis based on the global digital crop maps developed by Monfreda et al. (Monfreda et al. 2008) The nine major crop zones selected in this study are Central Europe and Russia (Eur_Rus) zone, East China zone, North America zone, South Australia zone, Southeast Latin America zone, West Latin America zone, Southwest Africa zone, Southwest Asia zone, and West Europe zone.Geographic locations of each major crop zone are given in Figure S1.
The statistics of main crop yields and harvested areas were downloaded from the FAO statistics division (http://faostat3.fao.org/), and the final yield of each major crop zone was computed through the harvested area weight-based average of yields in its corresponding inclusive countries.The inclusive countries of each major crop zone were listed in Table S1.

Method
The GIMMS dataset was first filtered by an automated compound smoother named RMMEH (Jin and Xu 2013) to efficiently reduce remaining noise and reconstruct high quality NDVI time-series dataset (Figure S2).A simple maximum value composite (MVC) method was applied to the filtered GIMMS data of each year to generate annual maximum NDVI data.
To detect the change trend of NDVI and climate variables (temperature, precipitation, and drought) over the entire study period , a least-square linear regression model (Piao et al. 2011) was applied as follows: where y represents annual NDVI or climate variable, t is year, a and b are the least-square fitted coefficients (a is the intercept and b is the trend slope), and ε is the residual bias.
A double subtraction algorithm (Peng et al. 2012) was employed to calculate the number of peaks in the filtered GIMMS temporal profiles of each year (Figure S3).The multiple cropping index (MCI) is used to indicate the planting frequency and intensity.
Climate change is regarded as one of major drivers of the crop productivity changes.For nine major crop zones, we applied the linear regression analysis between NDVI and climate variables to detect whether a positive or negative correlation exists.A p value < 0.05 was considered significant in this study.
Because of the regular latitude/longitude grid used in this study, the area of grid cells varies with the latitude.The global mean NDVI value or climate variables are therefore calculated through a grid-based weighted average as follow: where T is the grid-based weighted average result, ai is the grid area, and ti is the corresponding NDVI or climate variable value of the ith pixel.

Global crop growth trend over the past three decades
The max value of annual time-series NDVI can be used to define the magnitude of crop productions.Higher NDVI is expected to be an indicator of better agricultural growth or productivity.Figure 1a shows the spatial distribution of global crop maximum NDVI trend over the past three decades .A significant increased trend of crop maximum NDVI is observed globally.Statistically, 82.97% of global cropland witnesses a greener trend while 17.03% of that shows a yellower trend.With a zoomed major crop zone of the East China (Figure 1b), we can clearly identify the spatial distribution of greener and yellower crop growth.We further investigate the crop growth trend in terms of nine major crop zones (Figure S1).All the crop pixels of each major crop zone were averaged to produce the mean maximum NDVI value for each year.Figure 2a-b provide a visual assessment of crop growth trend of nine major crop zones from 1982 to 2012, based on the GIMMS dataset.Statistically, a significant upward trend of crop growth over the past three decades was observed for all global major crop zones, i.e., the Central Eur_Rus zone (R 2 = 0.68, p < 0.001), East China zone (R 2 = 0.24, p = 0.005), North America zone (R 2 = 0.58, p < 0.001), South Australia zone (R 2 = 0.29, p = 0.002), Southeast Latin America zone (R 2 = 0.73, p < 0.001), Southwest Africa zone (R 2 = 0.59, p < 0.001), Southwest Asia zone (R 2 = 0.62, p < 0.001), West Europe zone (R 2 = 0.64, p < 0.001), and West Latin America zone (R 2 = 0.66, p < 0.001).
Figure 1.Global crop maximum NDVI trend during the entire study period .(a) is the spatial distribution of global crop NDVI trend, and (b) is the corresponding zoomed sub-region (taking the complete geographical outline of East China Zone for an example).The green denotes the increased crop growth trend while the yellow represents the decreased crop growth trend.The darker green/yellow is expected to be a higher level of corresponding increased/decreased crop growth trend.
Coincident yield statistics of seven main crops from the Food and Agriculture Organization of the United Nations (FAO) were also integrated for an inter-comparison regarding each major crop zone.Figure 3 shows the annual yield change trend of main crops (i.e., barley, maize, rice, sorghum, wheat, and total cereals) during the entire study period 1982-2012.Visual inspections reveal that almost all crops witness an obvious upward yield trend in the nine crop zones (Figure 3a-i), and the corresponding statistical report in Table.S2 verifies that except the barley in the West Europe zone, the rice in the Southwest Africa zone, the sorghum in the Central Eur_Rus and North America zones, the Soybeans in the West Latin America and West Europe zones, the wheat in the South Australia and Southwest Africa zones, and the total cereals in the South Australia zone, the remaining crops' yields in all crop zones all achieve significant increasing trend (p < 0.05).Further analyses are performed on annual mean NDVI (Figure S4a) and annual maximum NDVI (Figure S4b) change trend along with different latitudes, and results show that both trends increase significantly for geographic latitude zones with 15° latitude intervals from 60° N to 45° S (p < 0.05).The corresponding statistic reports are provided in Table S3.

Spatiotemporal distribution change of multiple cropping system
Multiple cropping index (MCI) defines the planting frequency and intensity of the crop in the same arable land in one year, it is a composited reflection of the natural resources, including the utilization ration of water, soil, photosynthesis, and economic requirement corresponds to agricultural cropping systems.Figure 4 shows the spatial distribution of multiple cropping system in four selected years (1982, 1992, 2002, and 2012) with a 10-year interval.It indicates that single cropping system is the dominant cropping system globally, and cropping intensity increases gradually from areas with a higher to lower latitude, and for example, MCI increases from one to three from northern to southern China (Piao et al. 2010).The multiple cropping system including double or triple cropping system is mainly distributed in Asia (e.g., China, India, Thailand, Laos, and Vietnam), Southwest Africa, Southeast Latin America (e.g., Argentina), and Mexico.Increased crop intensity was observed during the past three decades in Figure 4a-d, and to be specific, the cropland area with double cropping system increased from 31.79% to 49.77% globally during 1982-2012.MCI change trend in terms of separate major crop zones in Figure 5 also reveals that almost all major crop zones witness an upward trend of average MCI change except a relatively stable MCI trend shown in three zones (East China in Figure 5b, West Latin America in Figure 5f, and Southwest Asia in Figure 5h).MCI in 1982MCI in , 1992MCI in , 2002 and 2012 with a 10year interval, respectively.The green represents the spatial distribution of single cropping system, the orange represents that of the twice cropping system, and the dark red represents that of the triple cropping system.The insets show the global percentage of different cropping systems.
Figure 6 illustrates the multiple cropping system percentage changes with respect to different latitude zones.It is quite interesting to identify that a significant upward trend of multiple cropping intensity is observed for N45°_60° (p = 0.015), N30°_45° (p = 0.050), and N0°_15° (p = 0.011) in the north hemisphere (Figure 6a), S0°_15° (p = 0.012), and S30°_45° (p = 0.040) in the south hemisphere (Figure 6b).The result indicates that the spatial distribution of multiple cropping system is expanded from zones with a lower latitude to a higher latitude, and zones with a higher latitude show possible increased multiple cropping intensity.MCI in 1982MCI in , 1992MCI in , 2002 and 2012 with a 10year interval, respectively.The green represents the spatial distribution of single cropping system, the orange represents that of the twice cropping system, and the dark red represents that of the triple cropping system.The insets show the global percentage of different cropping systems.

Climatic drivers of the crop NDVI change trends
Climate change is regarded as one of the major drivers of the vegetation or crop growth changes (Piao et al. 2011;Zhou et al. 2001).In this study, we partially stand on the side of climatic driven analysis of crop growth changes.The spatial pattern of climate variables' change trend during the period 1982-2012 was provided in Figure S5, including temperature, precipitation, and Palmer Drought Severity Index (PDSI).Statistic results show that global warming is quite significantly observed in Figure S5a, and to be specific, 91.59% of the terrestrial surface witnesses an upward trend of temperature globally.76.13% of the terrestrial surface is encountering increased precipitation while 23.87% of that is undergoing decreased precipitation (Figure S5b).47.65 % of the terrestrial surface shows a more serious drought (Figure S5c).The statistical report of climate change trend in terms of individual major crop zone is summarized in Table S4.
Using the linear regression analysis between NDVI and climate variables zone by zone in Table S5, we infer that crop NDVI is significantly positive correlated to the temperature for the Central Eur_Rus zone (p = 0.043), Southeast Latin America zone (p < 0.001), Southwest Africa zone (p = 0.017), Southwest Asia zone (p = 0.001), and West Europe zone (p = 0.012).The remaining zones including the East China zone, North America zone, South Australia zone, and West Latin America zone do not show significant positive correlation between NDVI and temperature.Nevertheless, as for the South Australia zone and West Latin America zone, it shows that NDVI is significantly positive correlated to the precipitation (p = 0.010, and p = 0.031, respectively).The statistically significant correlation between NDVI and precipitation is also observed for Southwest Africa zone (p = 0.010), and West Europe zone (p = 0.046).However, insignificant relationship between the crop NDVI and PDSI for any specific major crop zone is shown in Table S5 (p > 0.05).

DISCUSSION
Due to the relative coarse spatial resolution (~8 km) of GIMMS data, land cover changes related to cropland (e.g., from others to cropland, or cropland to others) under sub-pixel resolution during the past three decades are neglected in this study.However, in spite of some inevitable cropland related landcover changes, the general spatial distribution of main croplands does not change largely over time.It is reasonable to utilize the AVHRR GLCF data that corresponds to the same spatial resolution of GIMMS data as the benchmarking map to mosaic the cropland boundary during the entire study period .
Globally, trends in the GIMMS NDVI3g series compare favourably with trends derived from Landsat (Beck et al. 2011) and the Moderate Resolution Imaging Spectrometer (MODIS) NDVI (Fensholt and Proud 2012).Here we further designed to perform the inter-comparison between GIMMS NDVI3g and MODIS NDVI during the overlapped period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) in Figure S6, aiming at the efficiency check of long-time GIMMS dataset in monitoring global crop growth.Figure S6a indicates the trend of GIMMS NDVI profile achieves an acceptable agreement with that of MODIS NDVI data, in spite of some biases in the absolute value level derived from the sensor and calibration differences (Fensholt and Proud 2012).When plotting the time-series MODIS (MODIS/Time) significant linear regression slope values against time-series GIMMS (GIMMS/Time) significant linear regression slope values for cropland areas at a global scale (Figure S6b), it shows clear that the majority of pixels have a slope value lager than 0, and the density scatter shows the densest proportion also falls into the 1-to-1 areas with a high correlation between MODIS and GIMMS NDVI change trend.
The human-related activities including the social-economic policy management, land cover changes, and urban sprawl are also one type of major drivers addressing the crop growth change.With a visual inspection of the spatial distribution of yellower crop growth in Figure 1, an interesting phenomenon can be identified that areas with a yellower crop growth are mainly distributed the developing countries such as China, Mexico, Thailand, and so on.For these countries, their terrestrial surfaces are experiencing an unprecedented change during the past three decades caused by land-use activities including converting natural landscapes for human use and changing management practices on human-dominated lands (Foley et al. 2005).Taking the East China zone in Figure 1b for an example, the major yellow areas cluster in the Yangtze River delta, Pearl River Delta, and the Central China.These places are rapidly developing in China during the past decades, which results in the land reclamation at the sacrifice of croplands.On the other side, we only analysed the climate variables including temperature, precipitation, and drought, and its corresponding impact on crop growth in this study.In fact, some of the crop growth gains may also result from "Green Revolution" technologies, including high-yielding cultivars, chemical fertilizers and pesticides, and mechanization and irrigation (Foley et al. 2005;Matson et al. 1997).Besides, carbon dioxide (CO2) concentration is also another potential factor impacting on crop growth.Most plants growing in atmospheric CO2 higher than ambient exhibit increased rates of photosynthesis.Thus, taking the impact of human activities, and more compound climate-derived factors into consideration at both global and regional scales deems to be open topics for our further research.

CONCLUSIONS
This study sought to use long-term GIMMS dataset from 1982 to 2012 to provide the global crop monitoring in terms of crop growth and cropping intensity.The results show that 82.97% of global cropland maximum NDVI witnesses an increased trend while 17.03% of that shows a decreased trend over the past three decades.The spatial distribution of multiple cropping system is observed to expand from lower latitude to higher latitude, and the increased cropping intensity is also witnessed globally.In terms of regional major crop zones, results show that all nine selected zones have an obvious upward trend of crop maximum NDVI (p < 0.001), and as for climatic drivers, the gradual temperature and precipitation changes have had a measurable impact on the crop growth trend.However, no significant correlation between crop NDVI and drought index for any major crop zone is identified.Figure S3.Calculation of the MCI using double subtraction algorithm.The original NDVI data (blue lines) with 24 points in a complete year were first stacked in a column, then a dislocation phase subtraction with one temporal phase was applied to original time series (e.g., t2-t1, t3-t2), this produced a new temporal profile (green lines named as S1).For these 23 points, the positive value was assigned with value of +1, whereas the negative was assigned with -1.Subsequently, a repeated dislocation subtraction with one temporal phase was applied to the previously derived profile, producing a third temporal sequence (orange lines named as S2) with 22 points and the values of the data were -2, 0 and 2. Finally, the -2 points (red arrows) were reported as the peaks of this annual NDVI profile, and the total number of -2 points was nominated as the multiple cropping index (MCI).(a) is the MCI calculation flowchart of single cropping system, and (b) is that of double cropping system.

Figure 2 .
Figure 2. Maximum NDVI trend of nine major crop zones during the entire study period (1982-2012) using GIMMS dataset.

Figure 3 .
Figure 3. Annual yield change trend of seven main crops in nine major crop zones.

Figure 4 .
Figure 4. Spatial distribution of multiple cropping system in selected years from 1982 to 2012.(a)-(d) indicate the spatial distribution ofMCI in 1982MCI in  , 1992MCI in  , 2002   and 2012 with a 10year interval, respectively.The green represents the spatial distribution of single cropping system, the orange represents that of the twice cropping system, and the dark red represents that of the triple cropping system.The insets show the global percentage of different cropping systems.

Figure 5 .
Figure 5. Spatial distribution of multiple cropping system in selected years from 1982 to 2012.(a)-(d) indicate the spatial distribution ofMCI in 1982MCI in  , 1992MCI in  , 2002   and 2012 with a 10year interval, respectively.The green represents the spatial distribution of single cropping system, the orange represents that of the twice cropping system, and the dark red represents that of the triple cropping system.The insets show the global percentage of different cropping systems.

Figure 6 .
Figure 6.Changes in the percentage of multiple cropping system (double or triple cropping zones) along geographic latitude zones with 15o latitude intervals from 60 o N to 45 o S over the past three decades.(a) Changes in percentage of multiple cropping system in the north hemisphere, and (b) the corresponding changes in the south hemisphere.In (a) and (b) the scatters indicate the upward trends of zones N45 o _60 o (p = 0.015), N30 o _45 o (p = 0.050), N0 o _15 o (p = 0.011), S0 o _15 o (p = 0.012), and S30 o _45 o (p = 0.040) are statistically significant, while the insignificant trends are shown in zones N15 o _30 o (p = 0.797) and S15 o _30 o (p = 0.255).Noted that limited cropland is distributed above S45 o , thus we only show three zones including S0 o _15 o , S15 o _30 o , and S30 o _45 o in the south hemisphere.

Figure S2 .
Figure S2.Schematic diagram showing the use of RMEEH to reduce the noise of time-series NDVI dataset.(a) one-year GIMMS data with a biweekly interval; and (b) one-year MODIS data with a monthly interval.The blue lines are original NDVI while the orange ones are the filtered NDVI using RMMEH.Results show that RMMEH efficiently reduces noise of NDVI time series while preserving accurate temporal profiles.

Figure S4 .
Figure S4.Annual mean NDVI (a) and annual maximum NDVI (b) changes along geographic latitude zones with 15o latitude intervals from 60o N to 45o S over the past three decades.

Figure S5 .
Figure S5.Spatial distribution of climate change over the past three decades.(a) Annual average precipitation change trend, (b) Annual average temperature change trend, and (c) Annual average PDSI change trend.The insets show the percentage distribution of corresponding change trends, and the marks "++" and "--" represent the total percentage of increase and decrease trend, respectively.

Figure S6 .
Figure S6.The inter-comparison of change trend between GIMMS NDVI3g and MODIS NDVI during the overlapped period (2000-2012).(a) Annual maximum and mean crop NDVI change trend derived from GIMMS and MODIS datasets; and (b) Density scatterplot (global) of regression slope values from linear regression trend analysis of annual maximum MODIS NDVI and time period 2000-2012 against that of annual GIMMS NDVI and time period 2000-2012.