The Temporal Analysis of Regional Cultivated Land Productivity with GPP Based on 2000–2018 MODIS Data

The spatiotemporal change characteristics of Cultivated Land Productivity (CLP) are imperative for ensuring regional food security, especially given recent global warming, social development and population growth. Based on the hypothesis that the Gross Primary Productivity (GPP) is a proxy of land productivity, the Moderate Resolution Imaging Spectroradiometer (MODIS) data with 500-m spatial resolution and 8-day temporal resolution was employed by the Vegetation Photosynthesis Model (VPM) to calculate GPP in Jilin Province, China. We explored the level of CLP using the GPP mean from 2000 to 2018, and analyzed the changing trend and amplitude of CLP in the whole study period using both Theil–Sen median trend analysis and the Mann–Kendall (MK) test, and forecasted the sustainability of CLP with the Hurst exponent. The trend result and the Hurst exponent were integrated to acquire the future direction of change. The results revealed that: (1) The CLP level was generally high in the southeast and low in the northwest in cultivated land in Jilin, China. The area with the lowest productivity, located in the northwest of Jilin, accounted for 15.56%. (2) The majority (84.77%) of the area showed an increasing trend in 2000–2018, which was larger than the area that was decreasing, which accounted for 3.97%. (3) The overall change amplitude was dominated by a slightly increasing trend, which accounted for 51.48%. (4) The area with sustainability accounted for 33.45% and was mainly distributed in the northwest of Jilin. The area with anti-sustainability accounted for 26.78% and was mainly distributed in the northwest and central Jilin. (5) The Hurst exponent result showed that uncertain variation of CLP is likely to occur in the future over the entire region, and the central region is prone to display degeneration. Therefore, the results of this study indicated that quality improvement policy could be implemented for the middle-to-low yield fields in northwest Jilin, and dynamic monitoring and protection measures could be implemented for the areas with uncertain future changes and decreasing sustainability.


Introduction
The CLP plays a critical and irreplaceable role in maintaining regional sustainable development and ensuring national food security [1][2][3]. With increasing food demand, it is urgent to explore CLP dynamics spatially and temporally at national and regional scales [4]. However, many studies have been conducted to research the spatial distribution of CLP [5,6], such as land productivity potential and deteriorating environmental conditions in Jilin. Over the past 50 years, the average temperature has risen by 0.35°C/decade and the average precipitation has dropped by 16.34 mm/decade in this region [33]. In the last decade years, the organic matter content (OMC) of black soil decreased from 6% to 2% in most of Jilin, which has caused some black soil farmland to lose production capacity [34]. Fertilizer is being applied 10 times more than 30 years ago and the utilization of pesticides has increased from 2.2 tons to 25,000 tons, which is the main reason for soil compaction and soil pollution [35]. It is necessary to determine the changing trend and variation amplitude of CLP in Jilin under the background of climate change and soil degradation. Further, the prediction of changes in the short term is the key to optimize the allocation of cultivated land resources, to maintain the stability of productivity and enhance the low-and medium-yield fields in Jilin.
This paper is based on the hypothesis that the GPP is a proxy of land productivity, and in this paper, we analyze the CLP level, the trend and amplitude of change in Jilin from 2000 to 2018 and the future direction of CLP. The GPP time series is estimated by the VPM model derived from MODIS data. The level of CLP was estimated using the GPP mean from 2000 to 2018. The trend and amplitude of GPP were quantified by both the Theil-Sen median trend analysis and the MK test. The sustainability characteristic of future change was assessed by the Hurst exponent. The MK trend and Hurst exponent were integrated to express the future direction of change. The objectives were to determine: (1) the temporal and spatial characteristics of the CLP level; (2) the trend and amplitude of CLP during 2000-2018; and (3) the sustainability and variation direction of GPP in the short-term future. Based on this study, the level of regional CLP improvement and cultivated construction achievements are promoted by understanding the temporal variation and the spatial distribution of CLP in the Jilin.

Study Area
The Jilin Province is located in the middle of northeast China, from 121 • 38 to 131 • 19 E and 40 • 50 to 46 • 19 N, with a total area of 18.74 × 10 4 km 2 , as shown in Figure 1. Jilin Province has a temperate continental monsoon climate with four distinct seasons, with rain and heat in the same season. The study area is divided into the eastern mountainous area and the western plain by the central Black Mountain. The cultivated land is about 5.54 million ha, accounting for 28.89% of the Jilin Province area, 70% of which is distributed in the central and western plains. The cultivated land in the Jilin is mainly dry land and paddy fields. According to the difference in the geographical environment, the study area was divided into the northwest ecologically fragile area, the central area with high quality and contiguous cultivated land and the southeast area located in the mountainous area with scattered distribution. The main crops of Jilin Province are C 4 -based maize (Zea mays L.) and C 3 -based paddy rice (Oryza sativa L.).

Data Source and Pre-Processing
The 8-day composite MOD09A1 dataset of MODIS was downloaded from the USGS website (www.usgs.gov/). The MOD09A1 is a level 3 product from the Terra platform with 500 m spatial resolution, where geometrical corrections and cloud contamination removal were done. The quality assessment band-flagged pixels of high quality were used for the analysis. A total of 380 images from 1 May 2000 to 1 October 2018 were applied for GPP estimation. The HDF format was transformed into GeoTIFF format and the sinusoidal projection was converted to the WGS84 projection with the MODIS Reprojection Tool (MRT), which can be acquired from the NASA Land Processes Distributed Active Archive Center. The cultivated land grids with 500 m spatial resolution were collected from the agricultural land classification database. The images were masked by the cultivated land grids using the Google Earth Engine (GEE) platform. The meteorological datasets consisting of 8-day mean temperature data were collected from the National Meteorological Center of China. The data were from 51 meteorological stations between 2000 and 2018. The observations were interpolated over 8 days with the same resolution as the MOD09A1 for obtaining spatiotemporal distributions of temperature to estimate the GPP. The flowchart of this paper is shown in Figure 2.
Sustainability 2020, 12, x FOR PEER REVIEW 4 of 16 days with the same resolution as the MOD09A1 for obtaining spatiotemporal distributions of temperature to estimate the GPP. The flowchart of this paper is shown in Figure 2.     days with the same resolution as the MOD09A1 for obtaining spatiotemporal distributions of temperature to estimate the GPP. The flowchart of this paper is shown in Figure 2.

Methods
To provide temporal characteristics of the CLP in the Jilin Province, the following investigations were conducted: 1.
The VPM model was used to estimate GPP for C 3 -based paddy rice and C 4 -based maize. Then, two types of comparable GPP were obtained using an anomaly analysis and the GPP time series was generated.

2.
To investigate the temporal and spatial characteristics of CLP over a long period, the level of CLP was analyzed by GPP multi-years mean. The trend and amplitude were quantified by both the Theil-Sen median trend analysis and the MK test. 3.
The sustainability characteristic of future change was assessed by the Hurst exponent, which was integrated with the trend result to express the future direction of change.

GPP Estimation with the VPM Model and the GPP Mean Calculation
The VPM model is a GPP estimation model driven by the RS indexes, such as the EVI and land surface water index (LSWI), which was developed by Xiao [36] and is based on conceptual partitioning of vegetation canopies into nonphotosynthetic vegetation (NPV) and photosynthetically active vegetation (PAV) or chlorophyll (CHL). The EVI and LSWI were calculated by the MOD09A1 bands that include blue (459-479 nm), red (620-670 nm), near-infrared (841-875 nm) and shortwave infrared (SWIR, 1628-1652 nm). The calculation formulas of the EVI and LSWI and the calculation procedures of the model are referred to in [23]. Two indexes were combined with other datasets as input parameters of the VPM model to calculate the GPP time series data (Figure 3.). The formula is as follows.
where GPP is the cultivated land GPP (Kg·C/m 2 ), ε 0 is the maximum light use efficiency, and T scalar , W scalar , P scalar are the scalars for the effects of temperature, water, and leaf phenology on light use efficiency of vegetation, respectively. PAR (MJ/m 2 ) is the photosynthetically active radiation in a period, FAPAR PAV represents the fraction of PAR absorbed by the vegetation canopy.

Methods
To provide temporal characteristics of the CLP in the Jilin Province, the following investigations were conducted: 1. The VPM model was used to estimate GPP for C3-based paddy rice and C4-based maize. Then, two types of comparable GPP were obtained using an anomaly analysis and the GPP time series was generated. 2. To investigate the temporal and spatial characteristics of CLP over a long period, the level of CLP was analyzed by GPP multi-years mean. The trend and amplitude were quantified by both the Theil-Sen median trend analysis and the MK test. 3. The sustainability characteristic of future change was assessed by the Hurst exponent, which was integrated with the trend result to express the future direction of change.

GPP Estimation with the VPM Model and the GPP Mean Calculation
The VPM model is a GPP estimation model driven by the RS indexes, such as the EVI and land surface water index (LSWI), which was developed by Xiao [36] and is based on conceptual partitioning of vegetation canopies into nonphotosynthetic vegetation (NPV) and photosynthetically active vegetation (PAV) or chlorophyll (CHL). The EVI and LSWI were calculated by the MOD09A1 bands that include blue (459-479 nm), red (620-670 nm), near-infrared (841-875 nm) and shortwave infrared (SWIR, 1628-1652 nm). The calculation formulas of the EVI and LSWI and the calculation procedures of the model are referred to in [23]. Two indexes were combined with other datasets as input parameters of the VPM model to calculate the GPP time series data (Figure 3.). The formula is as follows.
where GPP is the cultivated land GPP (Kg·C/m 2 ), ε is the maximum light use efficiency, and Tscalar, Wscalar, Pscalar are the scalars for the effects of temperature, water, and leaf phenology on light use efficiency of vegetation, respectively. PAR (MJ/m 2 ) is the photosynthetically active radiation in a period, FAPARPAV represents the fraction of PAR absorbed by the vegetation canopy.  To obtain the spatial distribution of the level of CLP during 2000-2018, the GPP mean of the period was calculated by the following formula:

Meteorological
where is the mean during 2000-2018, and n is the number of years. To obtain the spatial distribution of the level of CLP during 2000-2018, the GPP mean of the period was calculated by the following formula:

Theil-Sen Median
where GPP is the mean during 2000-2018, and n is the number of years.

Theil-Sen Median
The Theil-Sen median was used to extract the trends of the GPP time series. Theil-Sen median trend analysis is a robust trend statistical method that uses the median slopes between all n(n − 1)/2 pair-wise combinations of the time series data, which contribute to reducing the influence of missing or abnormal data and is particularly effective for the estimation of the trend in a time series [37]. The slope of the Theil-Sen median represents an increase or decrease with the GPP series data between 2000 and 2018 on a pixel scale. The calculation of the slope was as follows: where i, j are the number of years in the study period, and GPP i and GPP j represent the GPP values at times i and j, respectively. An increasing trend is inferred from Slope GPP > 0, and a decreasing trend is inferred from Slope GPP < 0. The significance of the trend was quantified using the MK test.

The MK Test
The MK test, which is a non-parametric statistical method that does not require compliance with a certain distribution and is not disrupted by other abnormal values, was used to quantify the significance of the trend. In the MK test, the null hypothesis (H 0 ) explains that the sequence (GPP 2000 , GPP 2001 , GPP 2018 ) consisted of a sample of n independent and uniformly distributed random variables. The alternative hypothesis (H 1 ) of the two-sided test was that the distribution of x j and x k is not identical for all values of j and k. S was calculated as follows: where S is the test statistic with a normal distribution. GPP k are ranked from k = 1, 2, . . . , n -1 and GPP j from j = k + 1, . . . , n. Each data point of GPP k is used as a reference point and is compared with all other data points of GPP j : where the range of Z is (−∞, +∞), and |z| > µ 1−α/2 is the significance level, while the sign of this time series shows significant changes on the level of α. In this study, the α was set as 0.05 to analyze the significance of the GPP trend from 2000 to 2018 at a confidence level of 0.05.

The Hurst Exponent
The Hurst exponent was proposed by British hydrologist Hurst (1951), and later revised by Mandelbrot and Wallis (1969), for determining the sustainability of the time series in the future. This method was used to assess the sustainability of vegetation growth by detecting whether the NDVI has long-term dependence [37,38]. In this paper, the sustainability of CLP in the short term was judged by Sustainability 2020, 12, 411 7 of 16 the quantitative description of GPP dependence using the Hurst exponent. The calculation steps are as follows: (1) To define the time series GPP(t), t = 1, 2...n.
(2) To calculate the mean of the GPP time series, (3) To calculate the accumulated deviation, (4) To acquire the level difference, (5) To acquire the standard deviation sequence, According to previous research [39], the value of H ranges from 0 to 1, and there are three types: When H = 0.5, it indicates that the GPP time series is a random sequence without sustainability, and it is inferred that the trend of the CLP in the future is unstable. When H > 0.5, it means that the GPP time series is a sequence with sustainability, the future trend is consistent with past changes, and the trend of the CLP in the future is sustainability. When H < 0.5, it means that the GPP time series is a sequence with anti-sustainability, the future trend is opposite to past changes, and the trend of the CLP in the future is anti-sustainability.

The Characteristics of the Temporal and Spatial Distribution of the CLP Level
The level of CLP was accessed by removing the inter-annual variation difference of GPP from 2000 to 2018 (shown in Figure 4a,b). According to the classification of the natural breakpoint (Table 1), the CLP level was classified into levels 1-5, from low to high. The whole was 0.633 Kg·C/m 2 , which means it was level 3, called the provincial level. The temporal trend of the annual GPP during 2000-2018 is shown in Figure 5a. The whole area fluctuated in the range of 0.57-0.74 Kg·C/m 2 and exhibited an upward trend with a slope of 0.0063 Kg·C/m 2 yr in the monitoring periods. The maximum was in 2017 with 0.57 Kg·C/m 2 and the minimum was in 2000 with 0.74 Kg·C/m 2 . Temporal periodicity was seen with a five-year cycle of increasing and then decreasing. More specifically, the temporal change at different levels is shown in Figure 5b, where it is evident that all slopes were higher than 0, which means that they all showed an increasing trend from 2000 to 2018. However, the slope rate of each level was different, and the fastest change rate was 0.0086 Kg·C/m 2 yr in level 5. The slowest change rate was with an annual growth rate of 0.00427 Kg·C/m 2 in level 1. Other levels fluctuated within this range during this period. The important conclusion was drawn that a higher level of productivity had a slower rate of increase and the level with the lower productivity had a faster rate of increase.       The spatial distribution of the CLP level is shown in Figures 4 and 5c. The spatial distribution characteristics of the northwest at low altitude and southeast at high altitude were displayed in the cultivated land in Jilin. The level 3 area was the most abundant and accounted for 38.4% (Table 1). It was distributed in the central region, such as Siping City and Changchun City. The level 2 area accounted for 6.37%, which was the lowest, and was concentrated in Baicheng City and Songyuan City in the northwest region. The level 1 area occupied 15.56% of the total area, and was mainly concentrated in the northwestern part of Jilin, such as Baicheng City and Songyuan City. The level 4 area occupied 30.84% of the total area, and was mainly concentrated in the southeastern part of Jilin, such as Jilin City, Liaoyuan City, Tonghua City and Yanbian City. The level 5 area accounted for 8.8% of the total area. This area was mainly found alongside the level 4 area, including southeast Jilin City, north Tonghua City, Baishan City and Yanbian City.

The Trend and Amplitude of Change of the GPP from 2000 to 2018
Theil-Sen median trend analysis and the MK test were used to detect the trend and amplitude for each pixel. The results were divided into three categories according to the actual situation of the S GPP . There was no area where the S GPP was 0. Where the S GPP was between −0.0005 and 0.0005 were classified as stable areas. Where the S GPP was greater than or equal to 0.0005 were classified as increasing areas. Areas with an S GPP less than −0.0005 were classified as decreasing areas. Moreover, the significance test results of the MK test were used to assess the amplitude of change. At the confidence level of 0.05, the area was classified as significant variations (Z > 1.96 or Z < −1.96) or slight variations (−1.96 ≤ Z ≤ 1.96). By integrating the classification results of the Theil-Sen median trend analysis and the MK test, it was comparable to the data of trend variations of the CLP and was classified into five categories.
Areas percentages of CLP trends and their distributions are shown in Table 2 and Figure 6, respectively. (1) The areas where the CLP showed a significantly increasing trend occupied 33.29% of the total area. They were mainly distributed in Baicheng City and Songyuan City in the northwest region, but were also scattered in the west Siping City and northwest edges of Changchun City in the central region. (2) The areas where the CLP showed a slight increasing trend occupied 51.48% of the total area. These areas were distributed throughout the region and mainly concentrated in Changchun City, central Siping City, scattered in northwest Baicheng City, central Songyuan City in the northwest region and central Liaoyuan City and north Tonghua City in the southeast region. (3) The areas where the CLP showed a stable trend accounted for 11.26% of the total area. These areas were mainly found alongside the areas with a slightly increasing trend, including northwest Yanbian City and southeast Jilin City. (4) The areas where the CLP showed a slightly decreasing and a significantly decreasing trend accounted for 2.51% and 1.46%, respectively. These areas were scattered in Yanbian City, Liaoyuan City, and Tonghua City in the southeast region (Figure 6c). The amplitudes of all types are shown in Figure 7. The rates of significant changes, slight changes and no change were 0.00739 Kg·C/m 2 ·yr, 0.00398 Kg·C/m 2 ·yr, and 0.00113 Kg·C/m 2 ·yr, respectively. City and southeast Jilin City. (4) The areas where the CLP showed a slightly decreasing and a significantly decreasing trend accounted for 2.51% and 1.46%, respectively. These areas were scattered in Yanbian City, Liaoyuan City, and Tonghua City in the southeast region (Figure 6c). The amplitudes of all types are shown in Figure 7. The rates of significant changes, slight changes and no change were 0.00739 Kg·C/m 2 ·yr, 0.00398 Kg·C/m 2 ·yr, and 0.00113 Kg·C/m 2 ·yr, respectively.

The Sustainability and Direction of Change of the CLP in the Future
As shown in Figure 8, the average Hurst exponent H value of GPP in Jilin Province was 0.52, which was greater than 0.5 and indicates that the CLP change in the future will be consistent with the past. Namely, it will be sustainable. As far as the spatial distribution of the partition was concerned, the area with sustainability accounted for 33.45% and was mainly distributed in Baicheng City and Songyuan City in the northwest, and Liaoyuan City in the southeast, of Jilin province. An H value less than 0.5 suggested that the future direction of change of the CLP will be opposite to the past, namely, anti-sustainability, which accounted for 26.78% and was mainly distributed in the northwest of Baicheng City, Siping City and Changchun City in the central Jilin. The H value equal to 0.5 suggested that the future direction of change of the CLP will be random, namely, stable, which accounted for 60.23% and was mainly distributed in the northwest of Changchun City, Liaoyuan City and Jilin City in Southeast Jilin. Overall, the stable and anti-sustainability areas showed signs of large fluctuations in CLP that need to be monitored.
To reveal the direction of change of CLP in the future, the Theil-Sen S value, which characterized the past trend, was superimposed by the H value and characterized the sustainability of future changes to generate the comprehensive information of the trend and sustainability. The results were summarized into six categories from 'sustainability increased' to 'uncertain' as shown in Table 3 and Figure 9. The area shown to be uncertain in the future was the largest in the Jilin Province and accounted for 39.77%, which was mainly distributed in the Changchun City and Siping City in the central part. The 'sustainability increased' and 'anti-sustainability decreased' areas, which indicate that CLP will be improved in the short term, accounted for 25.61% and were mainly located in Baicheng City and Songyuan City. The 'sustainability decreasing' and 'anti-sustainability increasing' regions, where the CLP deteriorate in the short term, accounted for 28.52% and were mainly located in the southeastern City and Songyuan City in the central part. The area that was shown to be stable in the future accounted for 6.1%, and was scattered in the north of Baicheng City.

The Sustainability and Direction of Change of the CLP in the Future
As shown in Figure 8, the average Hurst exponent H value of GPP in Jilin Province was 0.52, which was greater than 0.5 and indicates that the CLP change in the future will be consistent with the past. Namely, it will be sustainable. As far as the spatial distribution of the partition was concerned, the area with sustainability accounted for 33.45% and was mainly distributed in Baicheng City and Songyuan City in the northwest, and Liaoyuan City in the southeast, of Jilin province. An H value less than 0.5 suggested that the future direction of change of the CLP will be opposite to the past, namely, anti-sustainability, which accounted for 26.78% and was mainly distributed in the northwest of Baicheng City, Siping City and Changchun City in the central Jilin. The H value equal to 0.5 suggested that the future direction of change of the CLP will be random, namely, stable, which accounted for 60.23% and was mainly distributed in the northwest of Changchun City, Liaoyuan City and Jilin City in Southeast Jilin. Overall, the stable and anti-sustainability areas showed signs of large fluctuations in CLP that need to be monitored. To reveal the direction of change of CLP in the future, the Theil-Sen S value, which characterized the past trend, was superimposed by the H value and characterized the sustainability of future changes to generate the comprehensive information of the trend and sustainability. The results were summarized into six categories from 'sustainability increased' to 'uncertain' as shown in Table 3 and Figure 9. The area shown to be uncertain in the future was the largest in the Jilin Province and accounted for 39.77%, which was mainly distributed in the Changchun City and Siping City in the

Discussion
In terms of different regions, the level of CLP, the amplitude of the trend and the future sustainability were combined with the geographical characteristic to explore the driving forces that are contributing to guiding the CLP improvement projects. In the northwestern region of Jilin Province, the CLP level is low at level 2 ( Figure 5b) and a significant increase was shown during 2000-2018, which will keep sustainability increasing in the future. The average annual precipitation in this area is 400-700 mm, which is less than the 700-900 mm in the southeast of Jilin. Most of the cultivated land was low-lying, which will degenerate to saline-alkali land and wind-blown land, where the ecological environment is fragile. It is also more vulnerable to natural disasters such as flooding, salinization and degradation, and the artificial irrigation facilities are poor in this area. These factors limited the regional CLP level. As shown in Figure 10, precipitation has greater regional differences compared to temperature in the northwest, than in the central part. Therefore, precipitation is more likely to affect regional productivity levels than temperature. In recent years, measures such as the construction of a cultivated land shelterbelt, high-standard cultivated land construction and watersaving irrigation projects have been implemented in the central and western regions to effectively control disasters and improve environmental conditions, which might be the reason for the significant increase in the region. The improvement of climatic conditions such as the increase of precipitation and the rise of temperature over the years is conducive to agricultural production, while the development of cultivated land regulation projects has gradually improved the water and

Discussion
In terms of different regions, the level of CLP, the amplitude of the trend and the future sustainability were combined with the geographical characteristic to explore the driving forces that are contributing to guiding the CLP improvement projects. In the northwestern region of Jilin Province, the CLP level is low at level 2 ( Figure 5b) and a significant increase was shown during 2000-2018, which will keep sustainability increasing in the future. The average annual precipitation in this area is 400-700 mm, which is less than the 700-900 mm in the southeast of Jilin. Most of the cultivated land was low-lying, which will degenerate to saline-alkali land and wind-blown land, where the ecological environment is fragile. It is also more vulnerable to natural disasters such as flooding, salinization and degradation, and the artificial irrigation facilities are poor in this area. These factors limited the regional CLP level. As shown in Figure 10, precipitation has greater regional differences compared to temperature in the northwest, than in the central part. Therefore, precipitation is more likely to affect regional productivity levels than temperature. In recent years, measures such as the construction of a cultivated land shelterbelt, high-standard cultivated land construction and water-saving irrigation projects have been implemented in the central and western regions to effectively control disasters and improve environmental conditions, which might be the reason for the significant increase in the region. The improvement of climatic conditions such as the increase of precipitation and the rise of temperature over the years is conducive to agricultural production, while the development of cultivated land regulation projects has gradually improved the water and fertilizer conservation capacity of cultivated land. At the same time, due to the continuous improvement in the technical production level, such as the cultivation of drought-resistant and high-yield seeds and complete cultivated land management technology, the CLP level in the northwest region presents a sustainable increasing trend in the future. In the central region of Jilin Province, the average annual productivity is level 3 and it showed a Sustainability 2020, 12, 411 13 of 16 slight increase during 2000-2018, whereas unsustainability decreased and the future was found to be uncertain. The precipitation in the central region is higher than that in the northwest (shown in Figure 10). The artificial irrigation facilities are sufficient to meet the regional needs. Therefore, the productivity level has been higher in the central region for about twenty years. With the increase in precipitation and temperature in recent years [40], productivity is gradually increasing slightly. The average annual productivity level in the southeastern part of Jilin is at level 4. There were a slight increase and stable trend for 2000-2018, and a decrease in sustainability in the southeast region. The cultivated land in the southeast is mainly scattered sloping land. The soil humus is rich and the OMC is high. The precipitation is the highest in the whole region, and so the GPP is the largest. As the area with the highest CLP level in the whole region, with the increase of precipitation and temperature in the future, it may exceed the suitable temperature for crop growth. At the same time, a broken field in the south makes it less resistant to disasters, which may lead to a decrease of CLP in the future.
Sustainability 2020, 12, x FOR PEER REVIEW 13 of 16 unsustainability decreased and the future was found to be uncertain. The precipitation in the central region is higher than that in the northwest (shown in Figure 10). The artificial irrigation facilities are sufficient to meet the regional needs. Therefore, the productivity level has been higher in the central region for about twenty years. With the increase in precipitation and temperature in recent years [40], productivity is gradually increasing slightly. The average annual productivity level in the southeastern part of Jilin is at level 4. There were a slight increase and stable trend for 2000-2018, and a decrease in sustainability in the southeast region. The cultivated land in the southeast is mainly scattered sloping land. The soil humus is rich and the OMC is high. The precipitation is the highest in the whole region, and so the GPP is the largest. As the area with the highest CLP level in the whole region, with the increase of precipitation and temperature in the future, it may exceed the suitable temperature for crop growth. At the same time, a broken field in the south makes it less resistant to disasters, which may lead to a decrease of CLP in the future. The 500-m MODIS data can give very accurate surface information. Moreover, the cultivated land in Jilin Province is large and concentrated, so there are less mixed pixels in the 500-m MODIS data. Therefore, it is reasonable to use MODIS to study cultivated land in Jilin. The GPP was selected as the indicator of CLP because, compared with NDVI, the estimated model is combined with water, temperature, and photosynthetic radiation effective coefficient indicators, which are the conditions inherent in the productivity of cultivated land. The method combining the Theil-Sen median slope analysis and the MK test is superior to linear regression analysis, as it can not only obtain trend information but also quantify the amplitude of the trend, and furthermore, it has an early warning effect on preventing the degradation of productivity and appropriately reducing cultivated land-use intensity. In addition, the time series data do not need to obey specific distributions, strong error avoidance, and significance level tests. The sustainability, shown by the H value, only indicates whether the future changes are consistent with the past or not. The direction of change is still unclear. Therefore, the future CLP trend was quantified by superimposing the trend analysis results and the Hurst exponent. The result of coupling these tests is that it reflects the impact of past trends on the future more scientifically.
Based on the GPP time series of Jilin Province from 2000 to 2018, the trend in spatial and temporal variation, the amplitude of change and the sustainability of future changes of CLP were studied. However, different crop types may exhibit different inter-annual variations. The main crop types in Jilin are rice, corn, soybean, sorghum, etc. Therefore, more detailed work is needed to explore the effects of different crop types on GPP temporal and spatial changes. Based on studying the spatial and temporal changes of GPP for many years, combined with multi-source data such as soil type, high-standard farmland construction data, and other social and climate data, we could explore the driving force analysis of the trend, and provide effective regulation for providing regional productivity degradation policy. While the changes in GPP are affected by climate change, human The 500-m MODIS data can give very accurate surface information. Moreover, the cultivated land in Jilin Province is large and concentrated, so there are less mixed pixels in the 500-m MODIS data. Therefore, it is reasonable to use MODIS to study cultivated land in Jilin. The GPP was selected as the indicator of CLP because, compared with NDVI, the estimated model is combined with water, temperature, and photosynthetic radiation effective coefficient indicators, which are the conditions inherent in the productivity of cultivated land. The method combining the Theil-Sen median slope analysis and the MK test is superior to linear regression analysis, as it can not only obtain trend information but also quantify the amplitude of the trend, and furthermore, it has an early warning effect on preventing the degradation of productivity and appropriately reducing cultivated land-use intensity. In addition, the time series data do not need to obey specific distributions, strong error avoidance, and significance level tests. The sustainability, shown by the H value, only indicates whether the future changes are consistent with the past or not. The direction of change is still unclear. Therefore, the future CLP trend was quantified by superimposing the trend analysis results and the Hurst exponent. The result of coupling these tests is that it reflects the impact of past trends on the future more scientifically.
Based on the GPP time series of Jilin Province from 2000 to 2018, the trend in spatial and temporal variation, the amplitude of change and the sustainability of future changes of CLP were studied. However, different crop types may exhibit different inter-annual variations. The main crop types in Jilin are rice, corn, soybean, sorghum, etc. Therefore, more detailed work is needed to explore the effects of different crop types on GPP temporal and spatial changes. Based on studying the spatial and temporal changes of GPP for many years, combined with multi-source data such as soil type, high-standard farmland construction data, and other social and climate data, we could explore the driving force analysis of the trend, and provide effective regulation for providing regional productivity degradation policy. While the changes in GPP are affected by climate change, human interference, and national or local land and agricultural policies, the drivers associated with the CLP should be explored further. These explorations will be helpful to fully understand the regional CLP, and local governments should adopt policies or projects that promote regional CLP.

Conclusions
In this paper, the GPP time series data from 2000 to 2018 were estimated by the VPM based on 500-m MODIS MOD09A1 data. The multi-year mean, the Theil-Sen Median trend analysis, the MK test and the Hurst exponent methods were used to explore the CLP level, the trend and the amplitude of change, and the future sustainability and future direction of change of the regional CLP. The following conclusions were drawn: (1) Overall, the CLP level in Jilin Province was characterized during 2000 and 2018 by a medium level GPP. The CLP of Jilin was generally high in the southeast and low in the northwest. Therefore, when the cultivated land protection and land improvement measures were implemented, the northwest region with low productivity should be given more consideration.
(2) Most of the regions showed an increasing trend throughout 2000-2018, which was larger than that of the decreasing trend. Most changes were slight. Higher productivity had a lower rate of increase and lower productivity had a higher rate of increase.
(3) The area shown to be uncertain in the future is located in Changchun City and Siping City in the central part. These areas should be monitored dynamically. The improving areas were located in Baicheng City and Songyuan City, and should be given protection measures. The deteriorating area was mainly located in the southeastern cities and the central part, and should undergo comprehensive quality improvement measures.
(4) We propose that the temporal analysis methods are suited for exploring the trend and forecasting the direction of change in the future. Because they are universal, the methods could also be applicable for a finer data source and for other areas in the world to explore future trends.

Funding:
This research was supported by the National Key R&D Program of China, grant number 2017YFF0206801-2.