Assessing the Long-Term Evolution of Abandoned Salinized Farmland via Temporal Remote Sensing Data

Salinization in arid or semiarid regions with water logging limits cropland yield, threatening food security. The highest level of farmland salinization, that is, abandoned salinized farmland, is a tradeoff between inadequate drainage facilities and sustainable farming. The evolution of abandoned salinized farmlands is closely related to the development of cropping systems. However, detecting abandoned salinized farmland using time-series remote sensing data has not been investigated well by previous studies. In this study, a novel approach was proposed to detect the dynamics of abandoned salinized farmland using time-series multispectral and thermal imagery. Thirty-two years of temporal Landsat imagery (from 1988 to 2019) was used to assess the evolution of salinization in Hetao, a two-thousand-year-old irrigation district in northern China. As intermediate variables of the proposed method, the crop-specific planting area was retrieved via its unique temporal vegetation index (VI) pattern, in which the shape-model-fitting technology and the K-means cluster algorithm were used. The desert area was stripped from the clustered non-vegetative area using its distinct features in the thermal band. Subsequently, the abandoned salinized farmland was distinguished from the urban area by the threshold-based saline index (SI). In addition, a regression model between electrical conductance (EC) and SI was established, and the spatial saline degree was evaluated by the SI map in uncropped and unfrozen seasons. The results show that the cropland has constantly been expanding in recent decades (from 4.7 × 105 ha to 7.1 × 105 ha), while the planting area of maize and sunflower has grown and the area of wheat has decreased. Significant desalinization progress was observed in Hetao, where both the area of salt-affected land (salt-free area increased approximately 4 × 105 ha) and the abandoned salinized farmland decreased (reduced from 0.45 × 105 ha to 0.19 × 105 ha). This could be mainly attributed to three reasons: the popularization of water-saving irrigation technology, the construction of artificial drainage facilities, and a shift in cropping patterns. The decrease in irrigation and the increase in drainage have deepened the groundwater table in Hetao, which weakens the salt collection capacity of the abandoned salinized farmland. The results demonstrate the promising possibility of reutilizing abandoned salinized farmland via a leaching campaign where the groundwater table is sufficiently deep to stop salinization.


Introduction
Approximately 19.5% of irrigated land worldwide is affected by salt, posing a threat to food security and economic development [1]. Secondary salinity, a direct result of human interaction with the land due to poor management practices, is a prominent land degradation process in agricultural land, especially in arid and semiarid areas with shallow groundwater depth [2,3]. An improper irrigation and fertilization schedule introduces an amount of dissolved salt into the waterlogged area, where salinization occurs due to the capillary ascent of groundwater with high salt content [4]. Salt-affected soil limits growth and reduces the yield of a wide range of crops [5,6]. As the degree of salinization increases Previous studies have mainly focused on limited types of land-use classification or evaluating the soil salinity degree from non-vegetative areas. Most of them failed to identify abandoned salinized farmland, which is crucial in gathering salt from surrounding cropping areas [3]. Wu et al. [10] detected abandoned salinized farmland using monotemporal imagery via supervision classification. However, mono-temporal data provided limited information for the target and may be influenced by the timing of detection. Therefore, a robust and accurate method that can detect abandoned salinized farmland and evaluate the soil salinity degree simultaneously from temporal remotely sensed data is needed to assess the salinization process over a larger area.
In this study, a method based on the concept of decision trees was developed to monitor the long-term dynamics of salinization in the Hetao Irrigation District, an irrigation system that is over two thousand years old in North China. The direction of the salinization process can be affected by both natural and human factors, while a shallower groundwater table with higher evaporation demands is prone to salinization. In contrast, deeper groundwater levels with sufficient infiltration water (e.g., rainfall and irrigation) show positive signs of desalinization [27][28][29]. Human-related activities may be the dominant factor in determining the direction of salinization for human-dominated land use [30]. However, sometimes the human impacts caused by socioeconomic factors can be spontaneous, making such impacts difficult to predict. For example, a large-scale shift of the cropping structure may disturb the previous balance of the water-salt system, thereby affecting the salinization process. Therefore, the objectives of this study were: (1) to assess the spatio-temporal patterns of the abandoned salinized farmland via remote sensing data backed by a robust and accurate method; and (2) to analyze the influence factors that affect the long-term dynamics of the salinization process in Hetao.

Study Area
Covering a total area of 1.12 × 10 6 ha, the Hetao Irrigation District (40.1 • -41.4 • N, 106.1 • -109.4 • E) is one of the largest irrigation systems in China (Figure 1). The study area is approximately 250 km in width from west (Wulanbuhe Desert) to east (Wuliangsuhai Lake) and 50 km in length from south (Yellow River) to north (Langshan Mountain). There are mainly five types of surface coverage in the Hetao Irrigation District recorded in the 2020 land cover products provided by GlobeLand30 (http://www.globallandcover.com/, accessed on 31 August 2021). The sand dune area is mainly distributed in the western part, accounting for approximately 7% of the total area. The area with water bodies is approximately 6%. Wuliangsuhai Lake, the largest lake in Hetao, is located in the eastern part of the study area and receives drainage water from the irrigation system. The vegetative area and residential area account for 77% and 8% of the total area, respectively. Spring wheat, maize, and sunflower are the major crops in Hetao, with areas of 4.9 × 10 4 ha, 3.0 × 10 5 ha, and 2.8 × 10 5 ha in 2020, respectively (http://tjj.bynr.gov.cn/, accessed on 31 August 2021). Due to the typical arid and semiarid climate, the Hetao Irrigation District has deficient annual average precipitation (ranging from 60 mm to 280 mm), which is mainly concentrated in July and August, while the potential evaporation reaches 2000-2400 mm per year. The large divide between crop evapotranspiration and actual rainfall provokes high dependence on irrigation water for agriculture.
irrigating the uncropped field after harvest (usually conducted from mid-Oct mid-November). An amount of water is applied in the autumn irrigation to leach salinity downward, which balances the annual soil salinity. Unfortunately, autum gation only leaches the salt-affected surface soil, and there may be an adverse effe overdose of irrigation is applied. The soil saline content in the over-irrigate bounces back in spring after the thawing of the soil water.

Methods
As it is challenging to detect abandoned salinized farmland directly from sensing imagery, a method that adopts a step-by-step strategy to exclude the oth types was used in this study. First, all pixels (a pixel is an n-dimensional vector wh the number of investigated satellite imageries) in the stacked time-series NDVI i were classified into 20 unlabeled classes via the K-means clustering algorithm In Hetao, the main irrigation method is flood irrigation. In response to the national water-saving irrigation policy, drip irrigation has been deployed in Hetao since 2015. However, most areas with drip irrigation are small demonstration areas. Approximately 4.79 billion m 3 of water per year (2020, Statistics Bureau of Bayannur City) is diverted into the Hetao area from the Yellow River. The enormous scale of irrigation has lifted the groundwater table of Hetao, which varies between 0.5 and 3.0 m annually. In the meantime, irrigation water transports salt into the groundwater, which makes the salinity of the groundwater increase. Capillary water with high saline content brings salt into the field, and thus, soil salinization occurs. An artificial drainage system was established to combat the soil salinization process in the Hetao Irrigation District. However, this system operates with low efficiency due to the gentle slope of the terrain in the Hetao area, by which merely thirty percent of the inlet salinity was drained out of the irrigation district. To balance the soil salinity in the field, local farmers deploy a typical irrigation pattern, termed autumn irrigation; that is, irrigating the uncropped field after harvest (usually conducted from mid-October to mid-November). An amount of water is applied in the autumn irrigation to leach the soil salinity downward, which balances the annual soil salinity. Unfortunately, autumn irrigation only leaches the salt-affected surface soil, and there may be an adverse effect if an overdose of irrigation is applied. The soil saline content in the over-irrigated land bounces back in spring after the thawing of the soil water.

Methods
As it is challenging to detect abandoned salinized farmland directly from remote sensing imagery, a method that adopts a step-by-step strategy to exclude the other land types was used in this study. First, all pixels (a pixel is an n-dimensional vector where n is the number of investigated satellite imageries) in the stacked time-series NDVI imagery were classified into 20 unlabeled classes via the K-means clustering algorithm. Then, representative final clusters corresponding to the main land types were manually selected to build a plant shape model. Second, the pixels of non-vegetative areas and water bodies were detected by merging the final clusters based on the predefined threshold; for the vegetative area, the shape-model-fitting method was used to classify each cluster into a specific plant type. Next, sand dunes and residential areas were detected using pixel-based thermal information and the salinity index, and the area of abandoned salinized land was detected after stripping sand dunes and residential areas from non-vegetative areas. Finally, the number of pixels for each land type was counted, and the area was equal to the pixel number multiplied by 30 m × 30 m, where 30 m is the spatial resolution of Landsat imagery. A brief workflow chart is shown in Figure 2, and details of the methodology are described in the following sections. representative final clusters corresponding to the main land types were manually selected to build a plant shape model. Second, the pixels of non-vegetative areas and water bodies were detected by merging the final clusters based on the predefined threshold; for the vegetative area, the shape-model-fitting method was used to classify each cluster into a specific plant type. Next, sand dunes and residential areas were detected using pixel-based thermal information and the salinity index, and the area of abandoned salinized land was detected after stripping sand dunes and residential areas from non-vegetative areas. Finally, the number of pixels for each land type was counted, and the area was equal to the pixel number multiplied by 30 × 30 m, where 30 m is the spatial resolution of Landsat imagery. A brief workflow chart is shown in Figure 2, and details of the methodology are described in the following sections.   [31]. Consequently, we deployed only six ETM+ imageries for 2012 when the TM instrument was retailed and the OLI was still not online. In addition, high spatial resolution imagery from Sentinel 2 was used to establish a regression model between the   [31]. Consequently, we deployed only six ETM+ imageries for 2012 when the TM instrument was retailed and the OLI was still not online. In addition, high spatial resolution imagery from Sentinel 2 was used to establish a regression model between the salinity index (SI) and measured soil electrical conductivity (EC). Detailed information on the satellite data used is shown in Table 1. Information about acquisition time and identifiers of used data is summarized in Table S1 (Supplementary Material). The thermal band of the downloaded Landsat level-1 product was radiometrically calibrated with the brightness temperature data using ENVI 5.3 software (Harris Geospatial, Boulder, CO, USA). For other bands, the surface reflectance data were retrieved from the level-2 product. The missing information in the ETM+ imageries was recovered using the interpolation toolbox of ArcGIS 10.3 software (Environmental Systems Research Institute, West Redlands, CA, USA). All thermal imageries were interpolated to 30 m to maintain consistency with VI data For the Sentinel 2 level-1C product (top-of-atmosphere reflectance), the atmospheric correction procedure was conducted by the Sen2Cor toolbox developed by ESA (https://step.esa.int/main/snap-supported-plugins/sen2cor/, accessed on 31 August 2021). The tiles of imageries were mosaicked and subsequently clipped to fit into the predrawn boundary of the study area. The cloud contamination problem is the main obstacle to deploying optical remote sensing products for monitoring Earth's surface. The simple cloud score algorithm was used to calculate the pixel-based cloud score as follows [32]: cloud score = clamp[min(term 1 , term 2 , term 3 , term 4 , term 5 ), 0, 1] (1) where the cloud_score represents the cloud likelihood. A larger value of cloud_score indicates a greater likelihood of a pixel being cloudy. ρ x is the top-of-atmosphere (TOA) reflectance, as shown in Table 1. T b denotes the brightness temperature. A threshold of 0.3 was used to mask clouds in this study.

Ground Data
Soil samples were collected from the surface at a depth of 0-5 cm. The geolocation for each sampling point was recorded by a GPS device (GPSmap 76, Garmin Company, Olathe, KS, USA) which has an acceptable accuracy (±1.45 m). A total of 133 soil samples were collected from 33 preselected observation sites from May 2018 to October 2019 ( Figure 1c). Ten grams of dried and crushed soil sample was weighed and 50 mL of deionized water was added, mixed well, and allowed to stand. After the solution was layered, a conductivity meter (DDSJ-308A, Leici, Yidian Co., Ltd., Shanghai, China) was used to measure the EC value of the supernatant. Then, repeated EC measurements were averaged for each observation site. Linking the soil sampling sites and Landsat imagery was inappropriate due to the high spatio-temporal heterogeneity in the sampling area. As a result, the measured EC values at each site were estimated from the corresponding SI in the high spatial resolution Sentinel 2 imagery collected on 22 April 2019, when the area was uncropped and unfrozen. The statistical crop areas of Bayannur city, including the planting area for spring wheat, corn, and sunflower, were available at the Statistics Bureau of Bayannur City (http://tjj.bynr.gov.cn/, accessed on 31 August 2021). Notably, most cropland in Bayannur City was in the Hetao Irrigation District. Therefore, the statistics of Bayannur city were used to validate the results.

Analyzing Temporal Remotely Sensed Imageries
Detecting and excluding vegetative areas was the first step in detecting abandoned salinized farmland. The traditional mono-temporal VI threshold-based method is sensitive to the selected threshold, and it may be influenced by the timing of the detection (e.g., flooded irrigation and the variability of the planting schedule may degrade its performance). Therefore, time-series VI was used to analyze the dynamics of the land cover from different years to pursue higher accuracy. Potential land types were classified by the K-means clustering method and the land types of interest were detected and merged via the shapemodel-fitting method. The irrigation strategy is determined by cropping structure, which may affect the evolution of abandoned salinized farmland. As a result, the regional cropping structure was additionally detected from clustering results via the shape-modelfitting method [33].

Land Cover Classification via K-Means Clustering
K-means clustering [34], an unsupervised method that accepts arbitrary shapes of input for an individual year, was used to preliminarily classify Earth's surface. The applicability of unsupervised methods against supervised methods in different situations is discussed in Section 5 This algorithm assumes that the samples have k cluster centers and optimizes the cluster centers by minimizing the within-cluster distances (objective function) as follows: Remote Sens. 2021, 13, 4057 8 of 22 where C i is the cluster group and µ i is the cluster center for Group i. The initial positions in the centroids were randomly selected from the samples. The variation in initial centroids introduced uncertainty into the final result. Therefore, the K-means algorithm was run 10 times with different centroid seeds in this study, and the best case of initial centroids was selected based on its rate of convergence (K-means clustering was deployed by the Python Sklearn package). After clustering, representative final clusters corresponding to main land types were manually selected to analyze their typical VI patterns ( Figure 3) and to build a plant shape model to distinguish specific vegetation types (Section 3.2.2). The temporal reflectance dynamics of distinct surfaces are diverse due to the significant differences in the spectral reflectance characteristics of different landscapes. For a planting area, the number of peaks of the VI trajectory was determined by the number of cropping cycles. There is only one growing season in the Hetao Irrigation District due to the extremely low temperature in the winter. As a result, the VI dynamics of a cropping area show a single peak at the maximum leaf area index (LAI) stage (for instance, spring wheat reaches the maximum greenness at the heading stage at approximately DOY 160, while summer maize and sunflower peak at approximately DOY 210). The bodies of water show a negative temporal NDVI throughout the year that makes them easy to identify. The water body area accounts for a small proportion of the total area of the study area, and the spatial variability of water quality is relatively large. As a result, the standard error of NDVI is relatively larger for the body of water. The non-vegetative land has an NDVI with low variability, and the forest shows a long duration in growing season. However, the final classification results were unstable for different random seeds, although the initial centroids were optimized. For this reason, a large k, which is far greater than the actual classes, was used to explore more potential land types. Then the threshold-based method was used to detect and merge the non-vegetative area and water body, and the "shape-model-fitting" method was used to detect and merge specific vegetation types from the final clusters. A higher k allows the algorithm to split the vegetative and non-vegetative areas into land-use types at a fine scale but at the cost of heavy computation. The value of k was determined by the trial-and-error method, and k = 20 was finally applied to make the uncertainty of the final detection of land type small. Then, the types of interest were identified and merged by predefined rules. Specifically, cluster centers with an annual average NDVI (average value of NDVI time series for each year) ranged from 0 to 0.2 were regarded as non-vegetative land, and water bodies were confirmed if the value was negative. For the vegetative area, a shape-model-fitting method is introduced in the next section to identify specific crops based on phenological metrics.  The time-series NDVI for each cluster center was fitted by the shape model (an asymmetric logistic curve in this study) to identify the cropping structure of crops (Equations (9) and (10)).

Detecting Cropping Structure Based on Phenological Metrics
The time-series NDVI for each cluster center was fitted by the shape model (an asymmetric logistic curve in this study) to identify the cropping structure of crops (Equations (9) and (10)).
where f (x) is the shape model with five fitting parameters: a, b, c, d, and k ( Figure 4a). The shape model represents the typical temporal VI pattern for a specific crop [33,35]. a refers to the base value of the NDVI. b denotes the difference between the maximal NDVI and a. c is the DOY when the NDVI reaches the peak. d and k control the wideness and the skewness of f (x). The two inflection points of f (x) are: where the left inflection point denotes the fastest growth stage and the right inflection point corresponds to the senescence stage. Here, we introduced the growth active period (GAP) calculated as follows to represent the growing season length for each plant.

Detecting Sand Dunes and Deserts Using Thermal Information
Different crop fields (such as wheat and corn and sunflower) and non-vegetative land (including sand dune area, bare soil area, residential area, and abandoned salinized farmland) in Hetao were successfully detected by the K-means clustering and shape-model-fitting methods. However, the methods failed to subdivide the non-vegetative land with a similar temporal VI pattern. Thermal information was used to mask the sand dune area (or desert) out of the retrieved non-vegetative land using their significant difference in brightness temperature (Figures 5 and 6). Evapotranspiration (ET) of sand dune areas is largely determined by rainfall, resulting in a high-level net radiation (Rn) partition to the sensible heat flux (H). In contrast, abandoned salinized farmland ( Figure 5) and residential areas ( Figure 6) show lower surface temperatures than sand dune areas due to the higher air-specific humidity and soil moisture. The optimal fitting parameters were retrieved by fitting f (x) on the discrete VI points via the least square error method using the "scipy.optimize" module in Python. The fitted shape models for the dominant land cover types are shown in Figure 4b. Spring wheat was planted ahead of other crops, and thus, parameter c was significantly smaller. The forest area shows a longer GAP than the crops, and an infinite GAP was observed for nonvegetative land cover types. However, distinguishing corn and sunflower is challenging due to the similar VI dynamics, which require data with a higher temporal resolution [13]. As a result, this study did not distinguish corn from sunflower but only classified them into one category containing plants with tall stems.

Detecting Sand Dunes and Deserts Using Thermal Information
Different crop fields (such as wheat and corn and sunflower) and non-vegetative land (including sand dune area, bare soil area, residential area, and abandoned salinized farmland) in Hetao were successfully detected by the K-means clustering and shape-modelfitting methods. However, the methods failed to subdivide the non-vegetative land with a similar temporal VI pattern. Thermal information was used to mask the sand dune area (or desert) out of the retrieved non-vegetative land using their significant difference in brightness temperature (Figures 5 and 6). Evapotranspiration (ET) of sand dune areas is largely determined by rainfall, resulting in a high-level net radiation (Rn) partition to the sensible heat flux (H). In contrast, abandoned salinized farmland ( Figure 5) and residential areas ( Figure 6) show lower surface temperatures than sand dune areas due to the higher air-specific humidity and soil moisture.   Figure 5. Example of the difference in brightness temperature between sand dune areas and abandoned salinized farmland. The top left subplot is the brightness temperature data (K) calculated from Landsat 8 OLI imagery (acquired on 10 July 2019). The top right is the corresponding composed RGB imagery. The red area of interest (AOI) denotes the sand dune area, and the green area is the abandoned salinized farmland. The lower subplot is the corresponding Google Earth map that shows the details of AOIs.  Otsu's method [36], which searches for the threshold by maximizing inter-class variance (between-class variance), was used to distinguish the sand dune area automatically from other non-vegetative surfaces. Specifically, the brightness temperature was calibrated from the Landsat thermal band (band 6 for TM and ETM+, and band 10 for OLI), and optimal threshold k* between two classes was iteratively calculated as follows: where σ 2 B is the between-class variance, and L denotes the maximum value of the investigated variable.

Detecting Abandoned Salinized Farmland
With a similar brightness temperature, the bare soil and residential areas were still mixed with the expected abandoned salinized farmland, although the sand dune area was masked out from the non-vegetative land. The abandoned salinity field showed a bright color in the composed RGB imagery (Figure 5), indicating that it reflected more radiation in the visual bands. The widely used SI calculated as follows was deployed to extract the abandoned salinity field [37,38].
where ρ blue and ρ red refer to the surface reflectance of the blue band and red band, respectively. The SI difference between abandoned salinized farmland and residential area is shown in Figure 7. Apparently, there was a higher-level SI in abandoned salinized farmland than in residential areas. Similar to the extraction of sand dune areas, the threshold-based method was used for the extraction. Finally, abandoned salinized farmland was detected, combining temporal VI data, thermal imagery, and SI data.
where ρblue and ρred refer to the surface reflectance of the blue band and red band, re tively. The SI difference between abandoned salinized farmland and residential a shown in Figure 7. Apparently, there was a higher-level SI in abandoned sali farmland than in residential areas. Similar to the extraction of sand dune area threshold-based method was used for the extraction. Finally, abandoned salinized land was detected, combining temporal VI data, thermal imagery, and SI data.

Performance Measures
In this study, a total of 150 validation points (including 22, 37, 23, 12, 15, 13, 1 12 points for wheat, corn/sunflower, other crops, forest, water body, sand dune, res tial area, and abandoned salinized farmland, respectively) were selected to valida accuracy of the land-use identification and abandoned salinized farmland identifi in 2019 with metrics of precision, recall, F1-scores, and accuracy. The validation p were labelled manually using very high spatial resolution imageries (satellite con tion from Planet Lab Inc. and Google Earth imageries). The spatial distribution o dation points is shown in Figure 1b. Precision (also called positive predictive value) fraction of relevant instances among the retrieved instances, while recall (also kno sensitivity) is the fraction of relevant instances that were retrieved. The F1-score c

Performance Measures
In this study, a total of 150 validation points (including 22, 37, 23, 12, 15, 13, 16, and 12 points for wheat, corn/sunflower, other crops, forest, water body, sand dune, residential area, and abandoned salinized farmland, respectively) were selected to validate the accuracy of the land-use identification and abandoned salinized farmland identification in 2019 with metrics of precision, recall, F1-scores, and accuracy. The validation points were labelled manually using very high spatial resolution imageries (satellite constellation from Planet Lab Inc. and Google Earth imageries). The spatial distribution of validation points is shown in Figure 1b. Precision (also called positive predictive value) is the fraction of relevant instances among the retrieved instances, while recall (also known as sensitivity) is the fraction of relevant instances that were retrieved. The F1-score can be understood as a combination of both precision and recall. It should be noted that accuracy takes into account only the total number of samples and thus may be a misleading metric for imbalanced datasets. In contrast, precision, recall, and F1-score consider the total number of samples for each class separately.

Evolution of Cropping Structure in Hetao
The variation in cropping structure during 1988-2019 is shown in Figure 8. The remote sensing results are in good agreement with the statistical data. The Hetao Irrigation District is included in the administrative region of Bayan Nur city. Consequently, the detection area is generally smaller than the statistical data (statistics include data outside the study area). In general, the results show that the cropland in Hetao has continuously expanded during the last thirty years (area of total cropland has increased from 4.7 × 10 5 ha to 7.1 × 10 5 ha). The spring wheat planting area initially accounted for nearly 50% of the crop planting area (1988). However, since 2000, the wheat planting area has shrunk year by year, and in 2019, it was only 6%. This is because an increasing number of local farmers have tended to plant economic crops instead (e.g., sunflower) to make higher profits. As a result, in 2000, the total area of maize and sunflower surpassed the area of wheat (2.0 × 10 5 ha), and in 2019, the area of maize and sunflower accounted for 70% of the total cropland.
year by year, and in 2019, it was only 6%. This is because an increasing number of lo farmers have tended to plant economic crops instead (e.g., sunflower) to make hig profits. As a result, in 2000, the total area of maize and sunflower surpassed the are wheat (2.0 × 10 5 ha), and in 2019, the area of maize and sunflower accounted for 70% the total cropland. The detected area for each crop was significantly underestimated in 2006 and 2 (residual error between statistics and estimates of the total cropland area up to 10.4 × ha and 10.1 × 10 5 ha, respectively). This is because the weather conditions of the stu area were cloudy in those years; therefore, there were fewer remote sensing data ava ble than in other years.  The detected area for each crop was significantly underestimated in 2006 and 2011 (residual error between statistics and estimates of the total cropland area up to 10.4 × 10 5 ha and 10.1 × 10 5 ha, respectively). This is because the weather conditions of the study area were cloudy in those years; therefore, there were fewer remote sensing data available than in other years.

Desalinization Progress in Hetao
A significant exponential relationship between the SI and measured soil EC values is shown in Figure 9. The EC map was estimated from the SI via the SI-EC relationship in Figure 9, and the salinity classes were determined by EC values using the widely used criteria as follows: (1) 0-4 dS/m: salt-free or slightly saline; (2) 4-8 dS/m: moderately saline; (3) 8-16 dS/m: highly saline; and (4) >16 dS/m: extremely saline [39]. The area was calculated by counting the number of pixels for each saline class ( Figure 10). Some research has indicated that the SI = √ ρ blue × ρ red shows the best results [38], while the results of that study were for the entire study area; different degrees of salinity were not studied separately. In the future, more samples should be collected to analyze different degrees of salinization separately. Only imagery collected in March and April was used to calculate the SI when the study area was crop-free and unfrozen. The missing data for specific years indicated that there was no cloud-free imagery available from March to April. The results show that the area of non-(or slightly) saline conditions continued to expand from 1988 to 2019. In contrast, the rate of decline for the moderately saline area has accelerated since 2000, while the highly saline area shrank sharply from 1993 to 2001. The area with extreme salinity decreased to a low level after 2000, indicating the construction of drainage facilities and promoting water-saving irrigation technology, which shows a positive effect on desalinization. specific years indicated that there was no cloud-free imagery available from March to April. The results show that the area of non-(or slightly) saline conditions continued to expand from 1988 to 2019. In contrast, the rate of decline for the moderately saline area has accelerated since 2000, while the highly saline area shrank sharply from 1993 to 2001. The area with extreme salinity decreased to a low level after 2000, indicating the construction of drainage facilities and promoting water-saving irrigation technology, which shows a positive effect on desalinization.  The spatial distribution of saline classes in Hetao from 1988 to 2019 is visualized in Figure 11. In 1988, most lands in Hetao suffered varying degrees of salinity, and the areas with extreme salinity were scattered throughout the study area. Due to the construction of drainage facilities, this situation was relieved in the next two decades. A pumping Figure 10. The dynamics of desalinization in Hetao: (a) salt-free/slightly saline and moderately saline; (b) highly saline and extremely saline. The area for each year was statistically analyzed using mono-temporal SI data from March or April when the study area was uncropped. The trends are depicted by moving average data with a window of 5.
The spatial distribution of saline classes in Hetao from 1988 to 2019 is visualized in Figure 11. In 1988, most lands in Hetao suffered varying degrees of salinity, and the areas with extreme salinity were scattered throughout the study area. Due to the construction of drainage facilities, this situation was relieved in the next two decades. A pumping station was set up at the outlet of the irrigation district (located at the easternmost boundary of Hetao) in 1975, and a field drainage system was built from 1988 to 1996, which drains water into Ulansuhai Nur Lake (east of Hetao). The facilities turned eastern Hetao salt-free (or slightly saline) (Figure 11b,c). However, the drainage system ran at low efficiency in western Hetao due to the flat terrain, which caused western Hetao to still be affected by severe salinity even in 2010 (Figure 11b,c). In the last decade, the changes in the cropping structure (planting more sunflower and corn instead of wheat) required less irrigation. This was because the deeper roots of sunflower and corn enhanced their ability to uptake water from the shallow groundwater table. In addition, the promotion of water-saving techniques reduced the inlet salt from the Yellow River. As a result, the saline situation was significantly mitigated in 2019, when most of western Hetao turned salt-free (or slightly saline) (Figure 11d).

Spatio-Temporal Pattern of Abandoned Salinized Farmland
The evolution of the detected area of abandoned salinized farmland from 1988 to 2019 is shown in Figure 12. The results show that both the area of abandoned salinized farmland and sand dunes have declined steadily in recent decades. The tense relationship between food demand and crop yield provoked the installation of irrigation systems (e.g., center-pivot sprinkler irrigation systems, Figure 6b) in the desert area. As a result, an increasing desert area was turned into an oasis in Hetao (area of desert decreased from 1.32 × 10 5 ha to 0.48 × 10 5 ha). In addition, the improved irrigation pattern and techniques, drainage system, and optimized cropping structure lowered the saline situation, allowing some of the abandoned salinized farmlands to be cropped again (the area of abandoned saline area was reduced from 0.45 × 10 5 ha to 0.19 × 10 5 ha). The spatial distribution of the abandoned salinized farmland in the last three decades is shown in Figure 13. Initially, the abandoned saline area was distributed relatively evenly across the study area. As increasingly highly saline land transformed into non-saline crop field, the abandoned salinized farmland vanished in both eastern and western Hetao. However, it can be seen clearly that there was a group of abandoned salinized farmland that has been in existence for 30 years in central Hetao. These wastelands are usually vast; therefore, it is difficult to leach the salt out of the fields. Moreover, the high potential evaporation in this wasteland made it a groundwater sink, gathering salt from the surrounding cropping area. The incoming salt ascended into the surface soil with the upward water flux, and it might have accumulated in the abandoned salinized farmland if effective leaching had failed to be deployed.

Spatio-Temporal Pattern of Abandoned Salinized Farmland
The evolution of the detected area of abandoned salinized farmland from 1988 to 2019 is shown in Figure 12. The results show that both the area of abandoned salinized farmland and sand dunes have declined steadily in recent decades. The tense relationship between food demand and crop yield provoked the installation of irrigation systems (e.g., center-pivot sprinkler irrigation systems, Figure 6b) in the desert area. As a result, an increasing desert area was turned into an oasis in Hetao (area of desert decreased from 1.32 × 10 5 ha to 0.48 × 10 5 ha). In addition, the improved irrigation pattern and techniques, drainage system, and optimized cropping structure lowered the saline situation, allowing some of the abandoned salinized farmlands to be cropped again (the area of abandoned saline area was reduced from 0.45 × 10 5 ha to 0.19 × 10 5 ha). The spatial distribution of the abandoned salinized farmland in the last three decades is shown in Figure 13. Initially, the abandoned saline area was distributed relatively evenly across the study area. As increasingly highly saline land transformed into non-saline crop field, the abandoned salinized farmland vanished in both eastern and western Hetao. However, it can be seen clearly that there was a group of abandoned salinized farmland that has been in existence for 30 years in central Hetao. These wastelands are usually vast; therefore, it is difficult to leach the salt out of the fields. Moreover, the high potential evaporation in this wasteland made it a groundwater sink, gathering salt from the surrounding cropping area. The incoming salt ascended into the surface soil with the upward water flux, and it might have accumulated in the abandoned salinized farmland if effective leaching had failed to be deployed.

Validation of Land-Use Classification and Abandoned Farmland Detection
A total of 150 verification points were used to evaluate the performance of the land-use classification and abandoned farmland detection. As seen in Table 2, all metrics reached at least 0.7 and above, within the acceptable range. The proposed method achieved outstanding performance for wheat and water body classification, with all metrics higher than 0.93. This is because the temporal VI patterns of wheat and water bodies are significantly different from those of other land uses, making them easier to classify. The confusion matrix in Figure 14 shows that pixels in residential areas were easily misclassified. This is because the components of the residential area are complex, which makes the pixels in them have a high degree of spatial variability. Overall, the proposed method yielded a feasible result with a kappa coefficient of 0.876.

Validation of Land-Use Classification and Abandoned Farmland Detection
A total of 150 verification points were used to evaluate the performance of the land-use classification and abandoned farmland detection. As seen in Table 2, all metrics reached at least 0.7 and above, within the acceptable range. The proposed method achieved outstanding performance for wheat and water body classification, with all metrics higher than 0.93. This is because the temporal VI patterns of wheat and water bodies are significantly different from those of other land uses, making them easier to classify. The confusion matrix in Figure 14 shows that pixels in residential areas were  easily misclassified. This is because the components of the residential area are complex, which makes the pixels in them have a high degree of spatial variability. Overall, the proposed method yielded a feasible result with a kappa coefficient of 0.876.

Discussion
Supervised machine learning algorithms have shown a solid capability to learn abstract features in nonlinear systems [40]. However, the input of these algorithms should

Discussion
Supervised machine learning algorithms have shown a solid capability to learn abstract features in nonlinear systems [40]. However, the input of these algorithms should be a uniformly formatted array. In other words, the shape of the input should be the same in the training and the testing phases [21]. Formatting the time-series input for remote sensing data with high temporal resolution (e.g., MODIS-based products) is applicable via cloud and cloud shadow compensation technology [41]. For high spatial resolution products (e.g., Landsat-based products), contamination from clouds heavily affects the consistency of formatted VI time series in different years (the number of available imageries varied depending on the climate condition of that year). Directly extending the short time-series data through time interpolation introduces large uncertainty in the testing phase, which hampers the performance of fine-scale surface classification by the supervised method. To test the applicability of the supervised algorithms for land cover classification using timeseries data with varying lengths (Landsat imageries), we examined one of the widely used algorithms, the random forest, to perform the task. The training samples were manually selected in specific years (years with more cloud-free data available) and tested in other years (results not shown). The trained model indeed performed well for years with more available time-series data. Unfortunately, it failed for years with short time series.
In contrast, unsupervised algorithms, such as K-means clustering, classify the samples by learning the distinct features from the samples themselves instead of the training set. Then, the clustering results are labelled automatically by prior rules. The advantages of K-means clustering make it deployable on any dataset with different lengths of time series. The results confirmed that this method showed stable performance and worked well in both cloudy years and less cloudy years.
A significant desalinization process was observed during the last thirty years, while the area of extremely saline conditions was reduced by approximately 50% and the area of highly saline conditions was reduced by approximately 70% (Figure 10). These results were deduced from the EC estimates of bare soil during the fallow season based on an established SI-EC relationship. However, the distribution of collected soil samples was imbalanced where the samples with extremely high EC were sparse. As a result, the prediction error of the established SI-EC relationship may be large in extremely saline areas. The only way to improve the robustness of the SI-EC relationship is to involve more available data. Therefore, more soil samples will be collected in the study area in future studies. In addition, Sentinel-2 imageries were used to build the EC estimation model because of the high spatial heterogeneity of the surface of the study area; therefore, remote sensing data with high resolution are needed to match the ground samples. However, the scale effect may reduce model performance when the SI-EC relationship is used in imagery with a lower resolution relationship (caused by the mixed pixels). Therefore, EC estimates may be more accurate on relatively uniform land with a smaller scale effect.
When a field is affected by high-level salinity, this reduces the yield of most crops and causes the field to be abandoned. Abandoned salinized farmland acts as an evaporative sink, resulting in a lower groundwater depth than the cropped area where constant irrigation increases the water table. Consequently, the salty water is drained from the cropped area into the abandoned salinized farmland, which mitigates the salinity in the cropping area. This is the concept of "dry drainage" [42]. The higher gradient of the water table in the dry drainage system transforms more salt into abandoned salinized farmland, which may increase its saline degree. The soil hardens in the abandoned salinized farmland when its saline degree is extremely high. The hardened soil layer reduces the capillary effect of soil water and decreases the efficiency of the dry drainage system. The tug-of-war between the efficiency of salt drainage and the evaporation capability of the abandoned salinized farmland creates a long-term salt balance in the dry drainage system [3]. In other words, the soil in the abandoned salinized farmland is difficult to improve due to the dry drainage effect, which has resulted in a thousand-year-old codependent relationship between the cropping system and the dry drainage system in Hetao.
The results show that the area of abandoned salinized farmlands has declined steadily in recent decades. The continuous decrease in the area of abandoned salinized farmland weakens the codependent relationship between the cropping system and the dry drainage system. There were three main reasons for the decline: (1) drainage facilities and watersaving irrigation techniques reduced the groundwater table [43]; (2) the inlet salt was reduced while the well-irrigation facilities were set up to pump groundwater for irrigation; (3) the introduced salt-tolerant crops that could grow in the salt-affected field not only lowered the area of the abandoned salinized farmland but also removed salt from the field [44]. The desalinization progress relaxed the tight connection inside the dry drainage system, where the deeper water table staggered the evaporation capability of the fallow area, and the shrinking abandoned saline area reduced the scale of the system. The soil salt content of wasteland decreased as the salt was leached by rainfall, while less salt was transferred into the area. The salt of the cropping area was directly leached into the deeper groundwater with irrigation and rainfall. This provided a moment to reutilize abandoned salinized farmland after artificial leaching campaigns. However, the risk of salinization bounce back still exists due to the rising groundwater table during rainy years. Therefore, a leaching project should be thoroughly designed and conducted step by step after evaluating the environmental conditions and the future climate.

Conclusions
In arid and semiarid water-logged areas, abandoned salinized farmland plays an essential role in acting as an evaporation sink gathering salt from the surrounding cropping area. However, traditional methods for land-use (land cover) classification fail to detect these wastelands. This study introduced a robust approach to detect abandoned saline land in the Hetao Irrigation District using 32-year time-series remotely sensed data. The proposed method is based on the idea of "elimination" to detect abandoned saline land step by step, which provides inspiration for monitoring the planting pattern of specific crops, or progress such as desertification or urbanization. At the same time, the dynamics of the cropping structure were monitored. The combination of the K-means cluster algorithm and the shape-model-fitting method successfully identified the specific crop, where the estimates agreed well with the statistics. The results show that the cropland in Hetao has continuously expanded during the last thirty years. An increasing number of farmers in Hetao have tended to plant high stem crops (maize and sunflower) instead of spring wheat from 1988 to 2019. The salinization degree was evaluated by the SI using an established exponential relationship between the SI and measured soil EC. The saline situation has been mitigated during recent decades, indicating that water-saving techniques, drainage facilities, and optimized cropping structures represent desalinization progress. Thermal information and the SI were used to distinguish the abandoned salinized farmland from other non-vegetative lands. The detected area of abandoned salinized farmland has continuously decreased, while many small wastelands have vanished in eastern and western Hetao. However, in central Hetao, a group of large-area abandoned salinized farmland has existed for 30 years. The strong effect of the dry drainage system has slowed the reutilization progress in these wastelands. However, a decrease in the groundwater table and salt-tolerant plants have created an opportunity to launch a saline field improvement project to reutilize the abandoned salinized farmland where the dry drainage system has already collapsed.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs13204057/s1, Table S1: Acquisition Time and Identifiers of Used Satellite Data.
Author Contributions: L.Z. performed data collection and processing; L.Z. and Q.Y. analyzed the data and interpreted the results; L.Z. and Q.Y. supported the writing, review and editing; J.W. and Q.Z. provided the suggestions and reviews. All authors have read and agreed to the published version of the manuscript.