Impacts of the world’s largest afforestation program (Three-North Afforestation Program) on desertification control in sandy land of China

ABSTRACT

In order to control the desertification, large-scale afforestation programs have been attempted worldwide. Among them, China initiated the world’s largest afforestation program, Three-North Afforestation Program (TNAP, 1978–2050), in which the afforestation in sandy land has been questioned during the first 40 years. In fact, the contribution of the TNAP to vegetation establishment and its effectiveness in desertification control still remain unclear, which limited the further construction of the program. To answer the questions, we detected the dynamics of vegetation distribution (forest, shrubland, and grassland) and desertification status (slight, moderate, severe, and extremely severe) during 1978–2017 in the sandy land (45.5 million ha), by visual interpretation of 5-period remote sensing images with validation based on 3,100 sample plots from field surveys and 15,175 sample plots from the National Forest Resource Inventory. Vegetation degradation was identified by analysis combining the trends of net primary productivity and precipitation use efficiency. By Geographical Detector model, the driving forces of vegetation degradation (climate change, human activities, vegetation type, and sandy land type) were ranked and the contributions of the influential factors (climate change, human activities, and vegetation dynamics) to desertification changes were estimated. The results showed that for the 40 years, vegetation coverage increased by 0.5%, with increasing 113.8% and 338.8% of forest and shrubland, but decreasing 9.0% of grassland. Desertification area had little change while the overall intensity decreased. The TNAP contributed to desertification dynamics by 10.3%, which is lower than expected. Vegetation type was the dominant factor of vegetation degradation in general. Forest is less suitable for afforestation in sandy land than shrubland and grassland because of its lower stand establishment rate, higher degradation rate, and less contribution to desertification control. Thus, adjusting vegetation type to match local conditions (e.g., use shrub-land, grassland, and native species) and improving the vegetation resistance (e.g., transform monoculture forests into mixed forests, and make proper proportion for forest, shrubland, and grassland) was suggested. Our study provided specific and feasible strategies for further planning and implementation of TNAP, and references for vegetation restoration of sandy lands worldwide.

KEYWORDS:

• sandy lands
• vegetation degradation
• desertification
• spatiotemporal changes
• driving forces

Introduction

Desertification is the most severe land degradation in drylands (including arid, semi-arid, and dry sub-humid regions) resulting from various factors, including climatic variations and human activities (UNCCD 1994). According to the report of The Millennium Ecosystem Assessment (Reid et al. 2005), desertification is one of the major environmental issues of the 21^st century in drylands because it would cause environmental and socio-economic-political problems at different scales, such as soil erosion, famine, and large-scale human migrations (Geist and Lambin 2004). Furthermore, the world’s drylands have significantly expanded over the past 60 years (Cherlet et al. 2018) due to the amplification of global warming (FAO 2019; Landman 2010). Since 1980, Asia has the world’s largest increases in dryland areas, of which, China alone accounts for about one-third of the increases exposed to aridity worldwide (Prăvălie et al. 2019). Vegetation degradation was the most direct distinguishing feature of land degradation and desertification (Bunning et al. 2016). Thus, to reverse the process of desertification, large-scale ecological restoration programs have been attempted worldwide (Gnacadja 2015; Chazdon 2008; Uzuner and Dengiz 2020).

China has a vast area of drylands in the Three North Regions (Northwest, North, and Northeast), in which sandified land is dominant, accounting for 172.1 million ha (17.9% of China’s territory area) (National Forestry and Grassland Administration 2014). Sandification is one form of desertification that characterized by the land surface being covered with sand and gravel, which often occurs in arid and semi-arid regions with fragile ecological ecosystems, or in areas around the deserts. The mobile sand driven by the wind encroach arable lands and pastures, ruin residents, and cause dust storms. This problem has been exacerbated by the increasing intensity of socio-economic human activities in the recent century (Zhang and Huisingh 2018). In order to control the expanding of desertification, the Chinese government launched the Three North Afforestation Program (TNAP), the world’s largest ecological engineering Program in 1978 and ongoing until 2050. The TNAP covered more than 42% of China’s land territory (over 4 million km²), including most of the arid and semi-arid zones of China, in which 96.5% of the sandy land of China (119.9 million ha) is located. During the past 40 years (1978–2017), the TNAP has made significant progresses in increasing local forest coverage and improving ecological conditions of the TNAP region (Wang 2014; Zhu and Zheng 2019). However, due to the water shortage, severe climatic conditions, and infertile soils, concerns have been raised about the effectiveness of afforestation in the TNAP (vegetation establishment), especially, in the sandy land of the TNAP region (Liu 2008; Cao et al. 2011; Zhang and Huisingh 2018). In fact, little is known about the contribution of the TNAP to vegetation establishment, and its effectiveness on desertification control in the first 40 years, which limits our understanding of how to manage the vegetation in sandy land, and the construction of the TNAP in the next three decades.

The estimation of vegetation establishment, i.e., the existence and status of vegetation, is hard to achieve for large ecological programs in drylands, due to the dynamic growth determined by inter-annual variations in meteorological factors, heterogeneity in geographical conditions, and difficulties in sampling for large areas in long time period (Verbesselt et al. 2010; He et al. 2019). Similarly, the process of desertification is the issue that should be studied on a regional or even national scale (Geist and Lambin 2004; Albalawi and Kumar 2013). The earth observation, i.e., satellite remote sensing, provided scientifically acceptable and powerful tools in monitoring vegetation dynamics and desertification process by recording changes in ecosystems on a long time and large space scale (Lanfredi et al. 2015; Smith et al. 2019). In landscape scale, the coverage area of the existing vegetation, i.e., the vegetation quantity, is the direct indicator in assessing the outcome of afforestation, which could be obtained from the interpretation and classification of vegetation from remote sensing images (Zheng et al. 2012; Zheng and Zhu 2017; Gerlein-Safdi et al. 2020). When assessing the vegetation status, i.e., the vegetation quality, the ecosystem structure, and function are the two aspects most concerned about (Smith et al. 2019). Based on remote sensing methodology, the vegetation index of net primary productivity (NPP) was considered to be more suitable for estimation of biomass (i.e., ecosystem structure) for sparse vegetation than other indices, such as normalized difference vegetation index (NDVI) (Anyamba and Tucker 2005; Sun et al. 2019). Since water is the limited condition in arid and semi-arid ecosystems, the ability to use available precipitation (PUE, i.e., precipitation use efficiency) is an apt feature in ecosystem functional assessment (Houerou 1984; Bai et al. 2008). Thus, as an integral method for estimating vegetation status (Li, Wei, and Zhou 2015), the combined analysis of NPP and PUE is especially suitable for identifying vegetation degradation in drylands.

As a process with distinct ground surface features, desertification could be identified through visual interpretation of remote sensing images with acceptable accuracy (Yan et al. 2009; Yan et al. 2015; Albalawi and Kumar 2013). However, the explanation for decades of desertification in spatial is a challenge. Over the years, previous studies have suggested that desertification is caused by several human factors (e.g., human activities, population dynamics, and agricultural policies) and natural factors (e.g., climatic change, soil type, and vegetation cover) (Du et al. 2017; Feng et al. 2015; Jiang et al. 2019). On that base, the contributions and the causal pattern of the factors are found to be specific to each locality (Geist and Lambin 2004; Yan et al. 2015; Sayl et al. 2021). Except for the irreducible complexity of desertification process, this finding also results from the differences in the methods of attribution (e.g., linear regression analysis and path analysis) and the difficulties in analyzing mixed types of data representing various factors (e.g., categorical data and continuous data). The model of Geographical Detector has been designed for quantifying the relative importance of several factors to response variables, free of the assumption of linear relationship of variables (Du et al. 2016; Wang et al. 2010). Based on the measurement of spatially stratified heterogeneity, the model can analyze mixed data types, including categorical data and continuous data (Cao, Ge, and Wang 2013; Wang, Zhang, and Fu 2016; Wang and Hu 2012). Having been successfully applied in the attribution analysis of vegetation dynamics at various scales (Zheng et al. 2021; Zhao et al. 2020; Zhang et al. 2021; Ran et al. 2019; Du et al. 2017), Geographical Detector is considered suitable for detecting the relative importance of vegetation to desertification dynamics among the other natural and human factors, and dominant drivers of vegetation degradation in sandy land of TNAP.

In this study, the study area was determined by extracting the sandy land from the dataset of Distribution of Desert (Sandy Land) in China (Wang et al. 2005) by excluding the desert (i.e., contiguous mobile sandy land and semi-mobile sandy land, referring to the high-resolution images provided by Google Earth) in the TNAP region, which is considered as the potential land for afforestation or vegetation restoration (45.5 million ha). Vegetation dynamics and desertification status were monitored and assessed by the interpretation of five periods of remote sensing images from 1978 to 2017 with the plot data of national inventory (National Forestry Administration 2000; National Forestry and Grassland Administration 2014, 2019, 2000; Xiao 2005; Jia 2009; Li et al. 2019) and field survey. Time-series data of vegetation biomass, climate, and human activities were used for studying the potential problems of existing vegetation and their driving forces. The specific objectives of this study were: (1) to display the vegetation dynamics both quantitatively and qualitatively; (2) to find the problems in present vegetation establishment and their driving forces; (3) to study the desertification processes and their relationship with the vegetation dynamics, climate changes and human activities, particularly the contribution of the TNAP to control the desertification processes. By achieving these objectives, specific and feasible strategies for further planning and implementation of TNAP are provided. References to vegetation restoration of sandy lands worldwide, especially in large scale and long term, are put forward.

Materials and methods

Study area

The study area was extracted by overlaying the distribution regions of sandy land and desertified land within the TNAP region (Figure 1). The sandy land was obtained from the Dataset of the Distribution of Desert (Sandy Land) in China (1:100,000) (Wang et al. 2005). The sandy land and deserts are mainly located in the northern provincial administrative units, including Xinjiang, Qinghai, Tibet, Ningxia, Inner Mongolia, Gansu, Shaanxi, Shanxi, Hebei, Beijing, Henan, Jilin, and Liaoning. The four main sandy lands (i.e., Hulunbuir Sandy Land, Otindag Sandy Land, Horqin Sandy Land, and Mu Us Sandy Land) of China are distributed in the east of Northern China. In the northwest, there are the eight main deserts of China (i.e., Taklimakan Desert, Gurbantünggüt Desert, Qaidam Basin Desert, Tengger Desert, Badain Jaran Desert, Kubuqi Desert, Ulan Buh Desert, Kumtag Desert) that consist of semi-mobile and mobile sandy land (Zhang and Huisingh 2018) (Figure 1a). The TNAP region is located in the Northeast, North, and Northwest of China, including the majority of sandy land and desertified land in China, in which the vegetation restoration is implemented for desertification control (Figure 1b).

Figure 1. Spatial distribution of study area. (a) shows the study area, consisted of sandy land and desertification area around the deserts in the TNAP; (b) shows all the sandy land and desert in the TNAP region (the four major sandy lands of Horqin Sandy Land, Otindag Sandy Land, Hulunbuir Sandy Land, and Mu Us Sandy Land are marked as I~Ⅳ, and the eight major deserts of Taklimakan Desert, Gurbantünggüt Desert, Kumtag Desert, Tengger Desert, Ulan Buh Desert, Kubuqi Desert, Qaidam Basin Desert, and Badain Jaran Desert are marked as Ⅴ~VII); and (c) shows the distribution of desertified land in the TNAP region.

The construction scope of the TNAP generally varied in different construction periods, and in the current period the TNAP has a total area of 435.8 million ha. To examine the interaction between vegetation restoration and desertification processes, the core area of deserts was excluded from the study area because the area is covered by stable desert. As a result, our study area covered 45.5 million ha, accounting for 4.7% of China’s land territory and 36.6% of total sandy land (124.2 million ha) of China (Figure 1c). According to the China’s Meteorological Background Dataset (Xu and Zhang 2017), the mean value of annual average temperature of the study area was 4.5°C (ranged between −6.6°C and 15.7°C), and the mean value of annual precipitation was 234 mm (ranged from less than 50 mm to around 470 mm). The original main natural vegetation types include desert steppe, typical steppe, meadow steppe, shrubland, and forest (Hou 2001; Fang et al. 2002; Su et al. 2020).

Data collection

The TNAP started in 1978. Remote sensing images of five time-points (1978, 1990, 2000, 2010, and 2017) were used to detect the quantity and spatial dispersion of vegetation and the status of desertification, representing the initial situation of vegetation and desertification before the TNAP, and the situation after 10, 20, 30, and 40 years’ implementation of the TNAP, respectively. Then the vegetation dynamics and desertification were analyzed in four periods (1978–1990, 1990–2000, 2000–2010, and 2010–2017). All the datasets used in this study are listed in Table 1.

Table 1. Overview of the datasets used in the study.

Data type	Dataset		Time period	Resolution	URL	Reference
Remote sensing images	Landsat MSS images	Landsat 3	1978	80 m	http://glovis.usgs.gov	(Loveland et al. 1996)
	Landsat TM images	Landsat 4, Landsat 5, Landsat 7 ETM+	1990, 2000, 2010	30 m
		Landsat 8 OLI	2017
	SPOT 5 images		2008	2.5 m	http://www.dsac.cn	(Latry and Rouge 1998)
	GF-1 images		2017	2 m, 16 m	http://www.cresda.com	(2015)
NPP	Global 5km8-day total and net primary productivity products		1981–2017	5 km	http://www.geodata.cn	(Cui et al. 2016; Tao et al. 2018; M. J. Wang et al. 2020)
Meteorological dataset	Monthly temperature and precipitation dataset for China		1981–2017	1 km	http://www.geodata.cn	(Peng et al. 2019)
	China’s Meteorological Background Dataset		2017 (date of publication)	1 km	http://www.resdc.cn	(Xu and Zhang 2017)
LUCC	China’s Multi-Period Land Use Land Cover Remote Sensing Monitoring Dataset (CNLUCC)		1980, 1990, 2000, 2010, 2015	1 km	http://www.resdc.cn	(Xu et al. 2018)
	Dataset of the distribution of desert (sandy land) distribution in China (1:100,000)		2005 (date of publication)	Shape file	http://westdc.westgis.ac.cn	(Y. Wang et al. 2005)
Socioeconomic dataset	China’s population spatial distribution dataset (1 kilometer grid)		1990, 1995, 2000, 2005, 2010, 2015	1 km	http://www.resdc.cn	(Xu 2017b)
	China’s GDP spatial distribution dataset (1 kilometer grid)		1995, 2000, 2005, 2010, 2015	1 km	http://www.resdc.cn	(Xu 2017a)

Remote sensing images

Five sets of remote-sensing images from 1978 to 2017 were selected, with Landsat TM images (30-m spatial resolution) acquired in 1990, 2000, 2010, and 2017, and Landsat MSS images (80-m spatial resolution) acquired in 1978, totaling 785 scenes (157 for each time point). The high-resolution images, i.e., SPOT 5 images in 2008 (2.5-m spatial resolution) and GF-1 images in 2017 (2-m spatial resolution), were used to improve the accuracy of vegetation monitoring and provide complementary information for the Landsat images, as in the previous studies of the TNAP (Yan et al. 2011; Zheng et al. 2012; Zheng and Zhu 2017), totaling 279 scenes and 1094 scenes, respectively. All the images were acquired between the end of June and mid-September, for the time period is the growing season of the local vegetation, and the images during the time period are apt for extracting the vegetation information (Zheng and Zhu 2017). To ensure the use of high-quality images (i.e., cloud free, good radiometric quality, low contents of aerosols, and high solar elevation angles), this study used the images of the previous year or the next year if the images of a specific year did not meet the requirements.

Satellite-based vegetation dataset

The spatial annual Net Primary Productivity (NPP) data from 1981 to 2017 were extracted from the Global 5 km 8-day Total and Net Primary Productivity Products (Cui et al. 2016; Yu et al. 2018; Wang et al. 2020). The dataset was produced by the light energy utilization model of MuSyQ-NPP based on 5-km resolution GLASS (Global Land Surface Satellite) LAI and FPAR products, and ERA-Interim meteorological data products. The model was improved by introducing a sky index, which could provide more accurate calculation of light use efficiency. The dataset was produced under the funding of the National Key R&D Program Project “Stereoscopic Observation and Inversion of Key Parameters of Global Ecosystem Carbon Cycle” and was verified by data from FLUXNET sites (http://www.resdc.cn).

Meteorological dataset

The meteorological data included moisture index, annual precipitation, and annual average temperature. The data of moisture index were obtained from the China’s Meteorological Background Dataset (Xu and Zhang 2017). The annual precipitation and annual average temperature from 1981 to 2017 were extracted from the 1 km monthly temperature and precipitation dataset (Peng et al. 2019). This dataset is downscaled in China through the Delta space downscaling scheme based on the global 0.5°Climate dataset provided by the Climatic Research Unit (https://crudata.uea.ac.uk/cru/data/hrg/) and the global high-resolution climate dataset released by the WorldClim (http://www.worldclim.org/), with verification from 496 independent meteorological observation stations throughout China.

Other dataset

The dataset of Distribution of Desert (Sandy Land) in China (1:100,000) (Wang et al. 2005) were used to extract the spatial range of sandy land. The local population and gross domestic product (GDP) were used as socioeconomic indicators, extracting from the China’s population and GDP spatial distribution dataset (1 km grid) (Xu 2017b, 2017a). The spatial distribution data of terrestrial ecosystem types, derived from the China’s Multi-Period Land Use Land Cover Remote Sensing Monitoring Dataset (CNLUCC) (1 km grid) (Xu et al. 2018), were used as a reference for the interpretation of remote sensing images. Data including National Forest Resource Inventory (the 5^th to 9^th) (National Forestry Administration 2000; National Forestry and Grassland Administration 2014, 2019, Xiao 2005; Jia 2009), Bulletin of National Desertification and Sandification Monitoring (the 1^st to 5^th) (National Forestry Administration 2000; National Forestry and Grassland Administration 2006, 2009, 2014; Li et al. 2019), national statistical data (construction/assessment reports of the TNAP and the Second National Land and Resources Survey Data), national permanent sample plot data (monitoring of ecosystem status and desertification, updated every five years), and field survey, were used to verify the interpretation results.

Field survey

The field surveys were carried out in the study area from June to September during 2007–2018 for investigation of vegetation and desertification status (Zheng et al. 2012; Zhu, Zheng, and Yan 2016; Zheng and Zhu 2017). 36033 Global Navigation Satellite System (GNSS) points were marked for interpretation and validation in the 375,165 km long survey route, with the minimum interval distance of 1 km. Land cover information, including vegetation type, vegetation coverage, and desertification intensity, were obtained by 3,100 sample plots, with each sample covering 900 m² of land area. The long-term monitoring data of the National Forest Resource Inventory (updated every five years) from 15,075 sample plots within the study area were used for the supplementary information of vegetation type, vegetation coverage, and vegetation health status, with each sample covering 667 m² of land area.

Classification of woody-grass vegetation and desertification

Image preprocessing

The 1:100,000 scaled digital topographic maps and other thematic maps (e.g., regional land-use, desertification, and vegetation maps) were used for geometric correction of the satellite images. All images were projected in Albers equal area conic projection, with geo-referenced and orthorectified by 70 to 80 uniformly distributed ground control points (GCPs) for each image. The root-mean-square error (RMSE) for the geometrical rectification using image-to-image matching method was at 1 to 2 pixels (i.e., <60 and <160 m for TM images and MSS images, respectively). All image data were processed by using an ERDAS 9.1 image processing system.

Visual interpretation and classification

The spatial distribution of vegetation and desertification at each time point were derived by the visual interactive interpretation method from the standard false-color images were composited by near-infrared band, red band, and green band (R, G, B) for Landsat MSS and TM remote sensing images (Yan et al. 2009; Zheng et al. 2012). Following the classification scheme developed in the previous studies of the TNAP (Zheng et al. 2012; Zheng and Zhu 2017), the vegetation of the study area was classified into five vegetation types at each time point (i.e., 1978, 1990, 2000, 2010, and 2017): coniferous forest, broadleaved forest, mixed conifer-broadleaved forest, shrubland, and grassland (Appendix 1). Among the vegetation types, the three arbor forest types and shrubland make up the planted forests in sandy land. Four intensity levels of desertification were identified from the false-color images, including slight desertification, moderate desertification, severe desertification, and extremely severe desertification (Yan et al. 2009; Zhu, Zheng, and Yan 2016), according to the classification system in Appendix 2. The land cover types (i.e., orchard, farmland, residential land, and waters) that were not considered in this study were included into a type named other land for statistical convenience. The vegetation types and desertification status were allocated manually according to the visual texture, color, and location of the patches in the composed images referring the natural vegetation maps and CNLUCC dataset (Zheng et al. 2012; Yan et al. 2009).

In the process of visual interpretation, the adopted mapping principles were listed as follows. First, the minimum map patch size was 6 pixels × 6 pixels, equivalent to 180 m × 180 m on the ground. Second, the positioning error was one pixel on the screen, equivalent to 2 m, 2.5 m, 30 m, and 80 m on the ground in the GF-1 images, the SPOT 5 images, the TM image, and the MSS images, respectively. Third, the classification results of vegetation type and desertification status were validated by the field investigation data, the National Forest Resource Inventory data, and other materials (e.g., regional land-use, desertification maps, and high-resolution maps from Google Earth), of which the overall accuracy was higher than 95% (Zhu, Zheng, and Yan 2016).

Calibration of vegetation quantity

Since the resolution of Landsat TM image is 30 m, the vegetation patch with areas less than 32,400 m² (by the mapping requirement at 1:100,000) was unable to be extracted. In order to improve the accuracy in estimation of vegetation quantity based on Landsat TM images, an area calibration method was put forward. Nine SPOT 5 images were randomly selected, of which the overlapped area with Landsat TM was split by plots with the size of 1 km × 1 km. In each plot, the vegetation area obtained from Landsat TM (i.e., medium-resolution) and by SPOT 5 (i.e., high-resolution) were extracted. Considering the variations in natural conditions across the study area, the plots were divided into two groups, i.e., annual average precipitation larger than 303 mm versus that lower than 303 mm (Zhu, Zheng, and Yan 2016). In each group, 80% of plots were used for the linear regression between the vegetation area obtained from Landsat TM and that by SPOT 5, and the other 20% were used for the verification. The accuracy of calibration was estimated as:

(1) $Pr=1-\left|\frac{1}{n}{\sum }_{i=1}^n\frac{x_e-x_0}{x_0}\right|\times 100\mathrm{\%}$(1)

where Pr (%) is the accuracy, xe is the calibrated value of vegetation area obtained from Landsat TM (m²), x₀ is the vegetation area obtained from SPOT 5 (m²), and n is the number of plots. The calibration relationship between vegetation area values extracted from SPOT 5 and Landsat TM used in this study is shown in Appendix 4.

Estimation of woody-grass vegetation quality

Calculation of precipitation use efficiency

Precipitation use efficiency (PUE) is a key indicator to assess vegetation status of dryland ecosystems, and has been widely used in studies on vegetation degradation and desertification (Houerou 1984; Hein and De Ridder 2006; Zhao et al. 2019). PUE is defined as the ratio between annual net primary production and annual precipitation, representing the ability of a vegetation system to convert precipitation water into vegetation (Prince, Brown De Colstoun, and Kravitz 1998). So, a decline in PUE means reduction in the capacity of the vegetation to transform water and nutrients to biomass, which indicates the vegetation degradation occurrence (Houerou 1984; Snyman 1998; Li, Wei, and Zhou 2015). In this study, the value of PUE for each time point was calculated as:

(2) $PUE=\frac{NPP}{P}$(2)

where PUE is the rain-use efficiency (g/m² mm) for each pixel, NPP (g/m²) is the average annual net primary productivity, and P (mm) is the annual precipitation.

Detection of the trends of NPP and PUE

The trend analysis was conducted for NPP and PUE in the four periods by acquiring the slope coefficient of linear regression at each pixel (1 km × 1 km). The slope coefficient is calculated as:

(3) $Slope=\frac{n{\sum }_{i=1}^{i=n}i\times x_i{\sum }_{i=1}^{i=\mathrm{n}}i\times {\sum }_{i=1}^{i=\mathrm{n}}{x}_i}{n{\sum }_{i=1}^{i=\mathrm{n}}{i}^2-{\left({\sum }_{i=1}^{i=\mathrm{n}}i\right)}^2}$(3)

where i is the number of the year in the period, and n is the range of the period, and $x_i$ represents the value of NPP or PUE at the ith year.

The significance of the slope was tested by the t-statistic. According to the outcome of the t-statistic, the pixels were classified into three categories: significantly increasing group (Slope >0, p < 0.25), significantly decreasing group (Slope <0, p < 0.25) and not significant group (p ≥ 0.25), referring to the method of Li et al. (Li, Wu, and Huang 2012). The trend of NPP and PUE was detected by four periods (i.e., 1978–1990, 1990–2000, 2000–2010, and 2010–2017).

Identification of vegetation quality types

Both NPP and PUE are considered to be effective factors to illustrate the vegetation degradation of arid and semi-arid ecosystem from different aspects (Prince, Brown De Colstoun, and Kravitz 1998; Hein and De Ridder 2006, 2006), and could be combined to improve the assessment of vegetation situation (Li, Wei, and Zhou 2015). In this study, degradation was detected by combining temporal change trends of NPP and PUE for the four periods, including 1978–1990, 1990–2000, 2000–2010, and 2010–2017. Based on the changes in both NPP and PUE, the existing vegetation was categorized into four quality types: severe degradation, potential degradation A, potential degradation B, and no degradation (Figure 2). For the two time points of Time 1 and Time 2, during the period between Time 1 and Time 2, significant decreases in both NPP and PUE indicated the vegetation degradation in both structure and function, which can be considered to be in a situation of severe degradation (Figure 2a). If NPP decreased while PUE increased or did not significantly change, the situation is considered as potential degradation A (Figure 2b). Under this scenario, total vegetation biomass decreased, but the vegetation self-recovery is possible because the function of the ecosystem did not decrease. The scenario that NPP increased or had not significant change but PUE decreased significantly is categorized into potential degradation B because the function of the ecosystem was actually deteriorating (Figure 2c). The situation that neither NPP nor PUE has a significant trend was divided into potential degradation B or no degradation. All the other situations, i.e., neither PUE nor NPP showing a significant decrease trend, are categorized into no degradation. Figure 2d showed one typical situation in which the structure and function of vegetation was improved. The outcome of the vegetation quality type identification was assessed by the National Forest Resource Inventory data, of which the overall consistency was around 80% (Appendix 5).

Figure 2. Four vegetation quality types identified based on the trends of NPP (net primary productivity) and PUE (precipitation use efficiency). The four vegetation quality types are (a) severe degradation, (b) potential degradation A, (c) potential degradation B, and (d) no degradation. The figure was adapted from the figure of Li, Wei, and Zhou (2015).

Statistical analysis

Spatiotemporal change of woody-grass vegetation

The land cover transition matrix (Table 2) was used to measure the spatiotemporal dynamics of vegetation during the four periods, i.e., 1978–1990, 1990–2000, 2000–2010, and 2010–2017. Time 1 and Time 2 refer to the starting and end points, respectively. $s_{ij}$ is the area of the land belonging to Type i in Time 1, but transfers into Type j in Time 2. If i equals to j, the land cover type remained unchanged from Time 1 to Time 2. $S_{n\ast }$ and $\hspace{0.17em}{S}_{\ast n}$ are the total area of the Type n in Time 1 and Time 2. During the period of Time 1 to Time 2, the decrease area of Type n equals to $S_{n\ast }-S_{nn}$, and the increase area of Type n is $S_{\ast n}-S_{nn}$.

Table 2. Land cover transition matrix during period of Time 1 - Time 2. This table was adapted from the table of Wang et al. (2019).

	Land cover type	Time 2				Total in Time 1
		Type 1	Type 2	…	Type N
Time 1	Type 1	$s_{11}$	$s_{12}$	…	$s_{1N}$	$S_{1\ast }={\sum }_{n=1}^N{s}_{1n}$
	Type 2	$s_{21}$	$s_{22}$	…	$s_{2N}$	$S_{2\ast }={\sum }_{n=1}^N{s}_{2n}$
	…	…	…	$s_{nn}$	…	…
	Type N	$s_{N1}$	$s_{N2}$	…	$s_{NN}$	$S_{N\ast }={\sum }_{n=1}^N{s}_{Nn}$
	Total in Time 2	$S_{\ast 1}={\sum }_{n=1}^N{s}_{n1}$	$S_{\ast 2}={\sum }_{n=1}^N{s}_{n2}$	…	$S_{\ast N}={\sum }_{n=1}^N{s}_{nN}$

Based on the transition matrix, net change area (NC), swap change area (SC), and total change area (TC) were calculated for the vegetation types. NC is the difference between the increase area and decrease area during the starting and end time points (Equation (4)).SC is the area of land that transferring into and out of one vegetation type (Equation (5)), which reflects the real change area during the process. TC is the sum of change area of one vegetation type during the period (Equation (6)). For vegetation type n, NC, SC, and TC were calculated as (Pontius, Shusas, and McEachern 2004; Wang et al. 2019):

(4) $N{C}_n=S_{\ast n}-S_{n\ast }$(4)

(5) $S{C}_n=2\times \mathrm{M}\mathrm{I}\mathrm{N}\left[\left(S_{\ast n}-S_{nn}\right),\left(S_{n\ast }-S_{nn}\right)\right]$(5)

(6) $T{C}_n=\mathrm{M}\mathrm{A}\mathrm{X}\left[\left(S_{\ast n}-S_{nn}\right),\left(S_{n\ast }-S_{nn}\right)\right]+\mathrm{M}\mathrm{I}\mathrm{N}\left[\left(S_{\ast n}-S_{nn}\right),\left(S_{n\ast }-S_{nn}\right)\right]$(6)

Geographical detector model

Geographical Detector Model (GeoDetector) was applied in the attribution analysis. GeoDetector is a method developed for exploring the relationships between landscape spatial patterns and their impact factors. It is a statistical tool for spatial variance analysis based on the measurement of spatial stratified heterogeneity (SSH), i.e., the phenomena that are within strata are more similar than are between strata (Wang et al. 2010; Wang and Hu 2012; Guo et al. 2022). Except that the assumption of linear relationship of variables is not required, GeoDetector is especially apt for analyzing mixed data types including categorical data and continuous data (Du et al. 2016; Liang and Yang 2016). The power of determinant (PD) of a certain factor is defined as:

(7) $PD=1-\frac{1}{N\ast {\sigma }^2}\sum _{i=1}^L{N}_i\ast {\sigma }_i^2$(7)

where $PD$ is the power determinant of the explanatory variable to target variable, $L$ is the strata number of the target variable stratified by the explanatory variable, $N$ and ${\sigma }^2$ are the total number of samples and the variance of target variable, respectively, $N_i$ and ${\sigma }_i^2$ denote the number of samples and the variance of target variable in the same stratum, respectively. $PD$ ranges from 0 to 1, representing that the explanatory power of an explanatory variable to the target variable goes from weak to strong. When $PD$ = 1, the target variable is completely dependent on the explanatory variable.

In this study, six factors were taken as potential explanatory variables into the driving factor analysis of vegetation degradation, including the trend of annual precipitation, the trend of annual average temperature, vegetation type, sandy land type, the change of local population and GDP, representing natural conditions and human activities. The six potential explanatory variables were calculated for the 10-year periods (i.e., 1990–2000, 2000–2010, and 2010–2017), in order to correspond to the results of vegetation quality types. In analyzing the driving forces for desertification process, area change for each vegetation type, the trend of annual precipitation and annual average temperature, the change of local population and GDP between 1978 and 2017 were included as explaining factors, with area change of each desertification intensity level as the target variable.

All the continuous factors (i.e., the trend of annual precipitation and annual average temperature, the change of local population and GDP, and the changes in area of vegetation and desertification) were discretized into nine levels by the quantile method in ArcGIS 10.4 for Desktop before putting into GeoDetector, because only discrete format explanatory variable data were applicable in the GeoDetector algorithm (Cao, Ge, and Wang 2013; Wang, Zhang, and Fu 2016; Liang and Yang 2016). The sample grid size for discretization in the driving factor analysis of vegetation degradation was set to be 1 km² to meet the resolution of majority of data. While 10 km² was taken as the sample grid size in analyzing the driving factors to desertification process, mainly considering that vegetation can affect the desertification process in a certain range around it. The GeoDetector software was downloaded freely from the website of http://www.geodetector.cn. The discretization and format transformation of the factor data were done by ArcGIS 10.4.

Results

Vegetation dynamics in sandy land during 1978 to 2017

The total area of vegetation in the study area increased from 30.9 million ha in 1978 to 31.2 million ha in 2017, with the coverage rate changing from 68.0% to 68.5% (Figure 3 and Appendix 6). Due to the afforestation programs (TNAP), forest, and shrubland increased by 113.8% (from 0.3 million ha to 0.7 million ha) and 338.8% (from 0.8 million ha to 3.3 million ha), respectively, during 1978 –2017. By contrast, grassland decreased by 9.0% (from 29.9 million ha to 27.2 million ha). Grassland was the dominant vegetation type in the study area at all time points, followed by shrubland and forest. Of the three forest types, broadleaved forests had the largest area, while conifer-broadleaved mixed forests had the smallest area.

Figure 3. Spatial distribution and statistic results of each vegetation type during 1978 to 2017. The spatial distribution of forest, shrubland and grassland in (a) 1978, (b) 1990, (c) 2000, (d) 2010, and (e) 2017 are shown. Each bar in (f) represents the area of the specific vegetation type of the year.

Vegetation type changes during the four periods of 1978–1990, 1990–2000, 2000–2010, and 2010–2017 were detected by the transition matrix (Appendix 6). For all vegetation types, the increase area was generally greater than the decrease area in the first three periods (1978–1990, 1990–2000, and 2000–2010), but contrary pattern was found during 2010–2017. Broadleaved forests had the largest positive net change area in the past 40 years among the three forest types, while mixed conifer-broadleaved forest was the only forest type that had a negative net change area. Shrub-land had the largest positive net change area during 2000–2010, but the largest negative net change area occurred during 2010–2017, resulting the positive net change area in the 40 years was much less than the total change area. Compared with the other three periods, the swap change area during 2010–2017 was the largest, indicating frequent and wide changes happened in vegetation during this period.

The distribution and proportion of the four vegetation quality types showed significant difference among the time periods (Figure 4 and Appendix 7). The overall degradation rate of forest and shrubland was 17.4% during 1978–1990, peaked at 40.5% during 2000–2010, and reduced to 25.3% during 2010–2017. For the forests, the degradation rate increased from 18.1% during 1978–1990 to 76.1% during 2000–2010, and reduced to 49.7% during 2010–2017. During 1990–2000, potential degradation A (degradation in structure) was the major degradation type, accounting for 23.7% of the total forest area (including 7.0 × 10⁴ ha of coniferous forests, 4.4 × 10⁴ ha of broadleaved forests, and 0.6 × 10⁴ ha of mixed conifer-broadleaved forests). In the next period, i.e., 2000–2010, potential degradation B (degradation in function) and severe degradation (degradation in both structure and function) were the two degradation types that happened most to the forests, accounting for 48.7% (including 8.3 × 10⁴ ha of coniferous forests, 16.6 × 10⁴ ha of broadleaved forests, and 1.2 × 10⁴ ha of mixed conifer-broadleaved forests), and 16.3% (including 5.6 × 10⁴ ha of coniferous forests, 2.8 × 10⁴ ha of broadleaved forests, and 0.1 × 10⁴ ha of mixed conifer-broadleaved forests) of total forest area, respectively (Figure 4 and Appendix 7).

Figure 4. Spatial distribution and statistic results of vegetation quality types for each vegetation type in four periods during 1978 to 2017. The spatial distribution of vegetation quality types for forest, shrubland and grassland during (a) 1978 to 1990, (b) 1990 to 2000, (c) 2000 to 2010, and (d) 2010 to 2017 are shown. Four vegetation quality types were identified, including severe degradation, potential degradation A, potential degradation B and no degradation. Each bar in (e) represents the proportion of area of the specific vegetation quality type of the period for each vegetation type.

The degradation rate of shrub-land fluctuated in the 40 years, changed from 17.1% during 1978–1990 to 24.6% during 1990–2000, peaked at 35.0% during 2000–2010, and reduced to 20.5% during 2010–2017 (Figure 4 and Appendix 7). Similar to those of forests, the dominant degradation type during 1990–2000 was potential degradation A, accounting for 10.8% of the total shrub-land area. During 2000–2010, the proportion of severe degradation shrubland did not increase significantly like the forests. Potential degradation A and potential degradation B were the two major degradation types during the time, accounting for 9.7% and 21.4% of the total shrub-land area, respectively. The results showed that for forest and shrubland, a high proportion of potential degradation B would happen following a high proportion of potential degradation A happened, indicating that the degradation in structure could result in the degradation in function in the following years (Figure 4 and Appendix 7).

Driving forces of vegetation degradation

The six potential driving forces affecting vegetation degradation were ranked according to their PD values (the power of determinant, Table 3). Vegetation type was the most important factor in the periods of 1990–2000 and 2000–2010, and the second most important in the period of 2010–2017. The highest rank of human activities factors, represented by the change of GDP (gross domestic product) and population, was 2^nd, 2^nd, and 3^rd during the three periods, respectively. Although the PD of temperature change was ranked as 1st during 2010–2017, the climate change factors were less influential than human activities in general (0.068 + 0.067 < 0.097 + 0.033 during 1990–2017). Sandy land type, i.e., the site conditions, had significant impacts on the vegetation degradation, although it was not ranked high.

Table 3. PD values and ranking of potential driving factors to vegetation degradation for period of 1990–2000, 2000–2010, and 2010–2017. PD values represent the percentage that each factor could explain the distribution of vegetation quality types in spatial in each time period.

Period		Climate change		Vegetation Type	Sandy Land Type	Human activities
		Precipitation Trend (mm/year)	Temperature Trend (℃/year)			Change of GDP (RMB/km^2·year)	Change of Population (person/km^2·year)
1990–2000	PD	0.010	0.014	0.062	0.010	0.038	0.018
	p-value	0.000	0.000	0.000	0.000	0.000	0.000
	Rank	6	4	1	5	2	3
2000–2010	PD	0.008	0.008	0.074	0.010	0.056	0.010
	p-value	0.000	0.000	0.000	0.000	0.000	0.000
	Rank	5	6	1	3	2	4
2010–2017	PD	0.001	0.045	0.013	0.004	0.003	0.005
	p-value	>0.1	0.000	0.000	0.000	<0.05	0.000
	Rank	-	1	2	4	5	3
Summation of 1990–2017	PD	0.018	0.067	0.149	0.024	0.097	0.033
	Rank	6	3	1	5	2	4

Desertification dynamics in sandy land during 1978 to 2017

Changes in desertification area and intensity

The total area of desertification was 24.9 million ha in 1978, peaked at 29.2 million ha in 2000 with an increase rate of 17.2% during 1978–2000, and reduced to 28.7 million ha in 2017 with a decrease rate of 1.5% during 2000–2017 (Figure 5 and Appendix 8). From 1978 to 2017, the area of all intensity levels of desertification generally increased with the time, except for extremely severe desertification, whose area decreased continuously after 2000. The proportion of slight desertification and moderate desertification to that of severe desertification and extremely severe desertification was 0.495: 0.505 in 1978, changed into 0.461: 0.539 in 2000, and became 0.551: 0.449 in 2017.

Figure 5. Spatial distribution and statistic results of desertification with intensity levels during 1978 to 2017. The spatial distribution of each desertification intensity levels in (a) 1978, (b) 1990, (c) 2000, (d) 2010, and (e) 2017 are shown. Each bar in (f) represents the desertification area of the specific intensity level of the year.

Driving forces of desertification dynamics

Seven factors were examined for their contribution to the desertification dynamics by PD values (Table 4). Change of vegetation area could explain 36.7% of the desertification dynamics during 1978–2017, while climate change and human activities explained by 12.5% and 25.3%, respectively. Of the vegetation factors, forest, shrubland and grassland explained 3.1%, 7.2%, and 26.4% of the desertification dynamics, respectively. Area changes at all intensity levels of desertification could be significantly explained by changes in forest, shrubland, and grassland areas. Grassland had the largest PD in all vegetation types for all intensity levels of desertification. For slight and moderate desertification, forest had similar explanation power than shrub-land (1.6% versus 1.5%). However, for severe and extremely severe desertification, the explanation power of forest was less than shrubland (2.7% versus 4.6%). Climate change and human activities had different explanation power in the desertification area changes for each intensity level (Table 4). Precipitation change explained the changes in slight and moderate desertification more than in the other two desertification intensity levels (4.5% versus 1.4%), while temperature change had the similar explanation power for all intensity levels of desertification (3.4% versus 3.3%). Human activities explained the area changes of desertification in slight and extremely severe levels (7.7% and 9.9%) more than the other two intensity levels (4.1% and 3.7%).

Table 4. PD values (%) of vegetation area change, climate change, and human activities to desertification dynamics during 1978 to 2017. PD values represent the percentage that each factor could explain the distribution of vegetation quality types in spatial in each time period.

Factor			Desertification Intensity Level				Summation PD of each factor
			Slight Desertification	Moderate Desertification	Severe Desertification	Extremely Severe Desertification
Vegetation area Change	Forest	PD	0.767	0.843	0.767	0.740	3.117	36.706
		p-value	0.000	0.000	0.000	0.000
	Shrubland	PD	0.816	1.837	2.255	2.333	7.241
		p-value	0.000	0.000	0.000	1.195
	Grassland	PD	4.072	6.410	10.400	5.465	26.348
		p-value	0.000	0.000	0.000	0.000
Climate change	Precipitation Trend	PD	3.154	1.313	0.541	0.852	5.860	12.477
		p-value	0.000	0.000	0.000	0.000
	Temperature Trend	PD	2.087	1.271	1.366	1.893	6.616
		p-value	0.000	0.000	0.000	0.000
Human Activities	Change of GDP	PD	2.081	1.537	1.425	4.608	9.652	25.288
		p-value	0.000	0.000	0.000	0.000
	Change of POP	PD	5.577	2.540	2.223	5.296	15.636
		p-value	0.000	0.000	0.000	0.000
Summation PD for desertification intensity level			18.554	15.752	18.978	21.187	74.470

Discussion

The impact of TNAP on vegetation restoration in sandy land of China during 1978 to 2017

Desertification is a combined consequence of multiple driving forces (Geist and Lambin 2004). The harsh natural conditions and the complex interactions between ecosystems and human activities in the sandy land of northern China present a challenge in the vegetation restoration through large-scale ecological programs (Zhang and Huisingh 2018). Our findings implied that vegetation restoration was not as effective as expected in the sandy land where the TNAP is implemented, since the vegetation coverage did not increase during the past 40 years (68.0% in 1978 and 68.5% in 2017). This is because the increase in forest and shrub-land by the TNAP is almost offset by the reduction of grassland. As reported by many previous studies, the low survival rate of the afforestation in the sandy land was believed to be one of the major causes for the unsatisfactory outcome of vegetation restoration (Wang et al. 2010; Wang et al. 2010; Cao et al. 2011). Stand establishment rate is another commonly used index in assessing the effectiveness of afforestation, which is the ratio of the actual increase area of forest or shrubland to the area of the newly planted forest or shrub-land that survived within 3–5 years. It was found that the establishment rates of the TNAP in sandy land was about 21.8% and 49.2% within 10 years, and decreased to 11.5% and 33.7% during the 40 years, respectively, for forest and shrubland. This indicates that only 1/5 of the survived forests became stands within 10 years of afforestation, and only 1/10 of the forests are still alive by 2017. In contrast, for the shrubs survived after afforestation, nearly half of them successfully became shrubland after 10 years, and more than 1/3 of the shrubland still exists by 2017. This suggests that the forest and shrubland suffered serious degradation consistently, but the shrub-land showed better adaptation than forest.

The increase of forest and shrubland area was limited not only due to the low establishment rate, but also due to the high degradation rate. It is obvious that the low quality or the degradation problem of the survived planted forests limited the functions of the ecological engineering projects (Li 2004; Zheng et al. 2012). It was found that the problem of vegetation degradation has always existed in sandy land during the construction of the TNAP. The degradation rate for forest and shrubland generally increased in the first 30 years, but decreased suddenly during 2010–2017. The possible reasons for this result might include: (1) the degraded forest and shrubland in 2000–2010 possibly died out during 2010–2017, as the decrease in area of both forest and shrub-land, during 2010–2017 was largest in all the periods. (2) There are likely some uncertainties in the assessment of degradation rates during 2010–2017, for the precipitation fluctuated more seriously in this period and the length of time for statistics was shorter than the other 3 periods. Therefore, the actual degradation rates during 2010–2017 may be larger than the calculated degradation rates in this study.

Among the factors that potentially caused the vegetation degradation in the study area, besides the mismatch between vegetation type and local site conditions, the socioeconomic factors were the dominant ones, except for the period when climate changed drastically (i.e., 2010–2017). According to Feng et al.’s study (2015), socioeconomic factors influenced more than environmental factors in NDVI decreasing and desertification increasing in China from 1983 to 2012. In previous studies conducted within the TNAP region, e.g., Altay Prefecture and Inner Mongolia, human disturbance was also proved to be the main factor responsible for vegetation degradation in the recent decades (Zhang et al. 2019; Sun et al. 2019). However, climate change became more influential than human activities in period with severe climate changes. With the dryer and warmer trend of climatic conditions in the northern China (Li, Chen, and Li 2019; Prăvălie et al. 2019), the meteorological factors would trigger more serious vegetation degradation in the future.

The contribution of TNAP in desertification changes in sandy land of China during 1978 to 2017

Through the construction of the TNAP, the desertification process during 1978–2017 in the sandy land was controlled to some extent, but the impact of the TNAP was very limited. Although the desertification area showed a decreasing trend after 2000, there was a net increase in the past 40 years. The overall intensity of desertification decreased, with slight and moderate desertification became dominant. The change trend of desertification area observed in this study was consistent with the conclusion from the 5^th Bulletin of National Desertification and Sandification Monitoring (National Forestry and Grassland Administration 2014).

The limited impact of the TNAP to desertification control was mainly due to the lower than preserved rate of forest and shrub-land (i.e., low survival rate and low stand establishment rate), and decrease in the area of grassland. The major determinants of the desertification dynamics in the study area during 1978–2017 was vegetation area change. The area changes of forest and shrubland explained 10.3% (3.1% and 7.2%, respectively) of the desertification dynamics, which was not as effective as expected. Grassland has more than double contribution to changes in desertification than forest and shrubland. Considering the fact that the area of grassland reduced by almost 9.0% in the study area, this may be one of the major reasons for the net increase in desertification area during 1978 to 2017. In addition, different vegetation types were found to have apparently different explanation of desertification area changes in each intensity level in this study. Of all the vegetation types, grassland worked well for dealing with all intensity levels of desertification. In controlling of slight and moderate desertification, forest worked as well as shrubland, while shrub-land worked much better than forest in controlling severe and extremely severe desertification. This is consistent with the findings in other studies that the establishment of shrub-land was an important step in the reversal of desertification processes (Maestre et al. 2009). The 40-year stand establishment rate of forest was only as 1/3 as that of shrubland and the forest had double degradation rate compared with the shrubland in the period with the most serious degradation, indicating that forest is less adaptive than shrub-land to the sandy land conditions. It was also found that the contribution of shrubland to desertification dynamics was more than twice of forest. These results suggest that shrub-land was more suitable than forest for large-scale afforestation in sandy land. These findings confirm that the future afforestation or vegetation restoration should fully consider the proportions of different vegetation types.

Human activities and climate change also weakened the impact of the TNAP in desertification control. This study reveals that the human activities were more crucial to desertification dynamics during the past 40 years than climate change. This result is consistent with the report that the socioeconomic factors were dominant in desertification development of China in recent decades (Li, Wu, and Huang 2012; Yan et al. 2015). The factors of human activities explained 25.3% of the desertification dynamics, took up 70.0% of the sum of the significant explanatory power by human activities and climate change (25.3% + 12.5%). In previous studies, this proportion was 65.8% to 68.6% (Zhou et al. 2014; Sun et al. 2019), which is similar to our results. Jiang et al. (Jiang et al. 2019) reported that the desertification process in grassland was more vulnerable to human activities than forest and sparse-forest-grass vegetation in Central Asia. In our study, it was clear that human activities not only resulted in the limited vegetation restoration, but also diminished the effect of the TNAP in desertification control in the sandy land. In most studies, the vegetation dynamics were used as a metric for indirectly identifying desertification rather than a driving force (Piao 2005; Verón et al. 2018; Jiang et al. 2019). Thus, our results provided a direct comparison between the effects of biological elements (i.e., vegetation dynamics) and non-biological elements (i.e., human activities and climate change) to desertification, which was an important supplement for understanding desertification. However, the intensity of human activities was difficult to describe adequately (Du et al. 2016). Thus the proxies used in this study, i.e., the change of population and GDP, may not be sufficient for detecting the full explanation power of human activities to desertification dynamics. According to Geist and Lambin’s study (Geist and Lambin 2004), desertification is driven by a limited set of variables but the pattern of synergies among causal factors has specific pathways for each region and period. The mixing of causal patterns caused by large space-time span of this study could be another reason for that the desertification dynamics was not fully significantly explained.

Suggestions and strategies

Here suggestions were proposed for the TNAP based on the findings in this study. For the future construction of the TNAP, the most important issue is to improve the survival rate and the stand establishment rate of afforestation, and minimize the degradation of surviving vegetation. According to our findings, the appropriate choice of vegetation type matching the local site conditions (particularly, water condition) could be the most effective way. With the originally harsh natural environment of the sandy land, the afforestation quantity and tree species should be determined according to scientific assessment of the site conditions, e.g., applying the annual precipitation, aridity index, or vegetation carrying capacity in planning (Zheng and Zhu 2017; Zhang et al. 2020). Besides the local species, shrub-land was demonstrated to be more adaptive to the natural conditions in sandy land than forest in this study. Thus, it is suggested to use shrubland, grassland, and native species to replace the forests in the areas of precipitation below 250 mm (Zheng and Zhu 2017). It is necessary to control human disturbances because the human activities were one of the major drivers for vegetation degradation and desertification. Grassland was found to be the most effective vegetation type in desertification control, the further disturbance of grassland should be urgently prohibited. According to the fact that the human disturbances, e.g., unsustainable grazing and cultivation practices, are related to the economic need of the local people, appropriate socio-economic strategies should be made based on “the interactions among the mountains, rivers, forests, farmland, lakes, grassland, sandy land, and the local people” (Yang 2004; Ge et al. 2015). Thus, rural poverty alleviation and the cost – benefit balance of the local community should be considered in the further construction of the TNAP (Cao, Suo, and Xia 2020; Suo and Cao 2021). Although climate change was less important than human activities in driving vegetation degradation and desertification, our study provides evidence that climate change could be a crucial factor of vegetation restoration in the future. With the stress from accelerated climate change (Huang et al. 2020), it is necessary to improve the resistance and resilience of the forest and shrubland under the conditions of climate change. For example, to establish the mixed conifer-broadleaved forests, because the mixed conifer-broadleaved forests were found to be able to keep lower degradation in all the four periods when other forest types suffered degradation. Similarly, the richness of plant species was also proved to be helpful in the functioning of restored vegetation (Li et al. 2021). According to the report of De Keersmaecker et al. (2015), the resistance and resilience of vegetation is also associated with the proportion of forest and non-tree vegetation (i.e., shrub-land and grassland). Thus, it is suggested that transforming monoculture forests into mixed forests in areas of precipitation more than 350 mm, and making proper proportion of forest, shrubland, and grassland in areas of precipitation more than 200 mm (Zheng and Zhu 2017).

There are uncertainties in the estimation of vegetation area in sandy land by remote sensing images in this study. First, errors were inevitable from the interpretation of remote sensing images. Even for the high-resolution images, e.g., SPOT 5 and GF-1, the vegetation patches with area less than 400 m² are difficult to be extracted (according to the pixel size of the images and 1:10,000 mapping requirement). Second, errors exist in the process of vegetation area calibration. Although the significant linear relationships were found in the vegetation area obtained from both Landsat TM and SPOT 5 (Appendix 3), there are still errors in the calibrated vegetation areas extracted from medium-resolution images compared with that from high-resolution images. Another limitation of this study is the method used for the identification of vegetation degradation by the combined trend analysis of NPP and PUE. In this study 10 years were used as the time period to detect the changing trend of NPP and PUE, because it is the minimum construction period of the TNAP. However, the 10-year period may not be the optimal time scale, as in other research the fluctuation period of NPP of vegetation in drylands was found to be 6 years, or even shorter (Jia, Liu, and Lin 2019). What’s more, the variations in NPP and PUE could be more complicated at the time scale of decades (Li, Wei, and Zhou 2015), which may not be able to be detected by linear regression. In future study, the trend analysis and degradation identification will be improved by combining the construction time point with the natural trend break point achieved by the Mann-Kendall method. In addition, to detect the spatial and temporal characteristics of NPP and PUE more precisely, multi time-scale analysis or multiple-regression will be applied. Besides, for research focused on large space-time scale like this study, the improvement of the quality of interpretation results is essential. In cloud platforms such as Google Earth engine, time-series composite remote sensing images could be managed, optimized, and compared, by which the seasonal effects and year-to-year change could be observed and captured better (Praticò et al. 2021). Thus, in future study, Google Earth engine will be applied for the improvement of vegetation classification.

Conclusions

This study revealed the vegetation dynamics and desertification process in the past 40 years in the sandy land with the construction of the TNAP. Our data showed that vegetation coverage increased by only 0.5%, with increasing 113.8% and 338.8% of forest and shrubland, but decreasing 9.0% of grassland during 1978–2017. Desertification had little change in area, but the proportion of slight-moderate intensity increased and outweighed that of the severe-extremely severe intensity. The TNAP contributed to desertification process by 10.3% in the 40 years, which is lower than expected. It was found that compared with shrubland, forest is less suitable for large-scale afforestation in sandy land, for it has relatively lower stand establishment rate, higher degradation rate, and less contribution to desertification reduction. In order to achieve the successful afforestation and more efficient desertification control, specific suggestions based on our research were provided as follows. (1) It is necessary to adjust the vegetation types to match the local conditions, e.g., using shrub-land, grassland, and native species to replace the forests in the areas of precipitation below 250 mm. (2) The measures to improve the vegetation resistance as well as the total impacts of vegetation on control of desertification should be implemented, e.g., transforming monoculture forests into mixed forests in areas of precipitation more than 350 mm, and making proper proportions for forest, shrubland, and grassland in areas of precipitation more than 200 mm. (3) The human disturbances should be controlled with appropriate socio-economic strategies by considering “the interactions among the mountains, rivers, forests, farmland, lakes, grassland, sandy land and the local people.” Our study provided a scientific support for planning and implementation of TNAP in the following 30 years and for the sustainable development of ecological restoration programs in sandy land worldwide. This study also provided an example for remote sensing research of sandy land vegetation on a large time-space scale, and further updates in methodologies combined with cloud platforms are suggested correspondingly.

Acknowledgement

This work was supported by the Consultation Project supported by Divisions of the Chinese Academy of Sciences [2019-ZW09-A-032]; and the National Nature Science Foundation of China [No. 31770758 and No. 32071837]. We acknowledge for the data support from “National Earth System Science Data Center, National Science & Technology Infrastructure of China (http://www.geodata.cn)”, “Resource and Environmental Science and Data Center, Institute of Geographic Sciences and Natural Resources Research (http://www.resdc.cn)” and “Cold and Arid Regions Sciences Data Center at Lanzhou (http://westdc.westgis.ac.cn)”.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request (http://www.iae.cas.cn/cslst/), except for the data from National Forest Resource Inventory and National Desertification and Sandification Monitoring.

Additional information

Funding

The work was supported by the National Nature Science Foundation of China [No. 31770758]; Consultation Project supported by Divisions of the Chinese Academy of Sciences [2019-ZW09-A-032].

Appendices

Appendix 1. Vegetation type classification system. This table was adapted from the tables of Zhu et al. (2016) and Zheng and Zhu (2017).

Land cover type		Description
Forest	Coniferous forest	Coniferous natural forests and plantation with crown density >30%, and height >2 m.
	Broadleaved forest	Broad-leaved natural forests and plantation with crown density >30%, and height >2 m.
	Mixed conifer-broadleaved forest	Mixed coniferous-broadleaved natural forests and plantation with crown density >30%, and height >2 m.
Shrub-land		Bush fallow and coppice with crown density >40% and height <2 m.
Grassland		Meadow steppe, real steppe and desert steppe with vegetation coverage >20%

Appendix 2. Desertification classification system and landscape features of each intensity level. This table was adapted from the table of Yan et al. (2009).

Desertification intensity level	Proportion of mobile sand or wind erosion area (%)	Description of landscape features
Slight desertification	<5	Small, sparse, scattered patches of mobile sand or wind-eroded areas appeared, with wind erosion pits on the windward side of sand dunes
Moderate desertification	5-25	Mobile sand sheets, coppice dunes, and eroded land scattered in the whole area, with obvious wind erosion slopes and falling sand slopes on sand dunes
Severe desertification	25-50	Semi-anchored dunes covered most of the area, with scattered mobile sand dunes or sand sheets. Eroded areas or sand sheets widely distributed
Extremely severe desertification	>50	Dense, mobile sand dunes or wind-eroded areas widely distributed

Appendix 3. Calibration relationship between vegetation area values extracted from SPOT 5 and Landsat TM and comparison of the vegetation area value before calibration and the value after calibration. The calibration relationship was built up separately for regions where annual average precipitation higher than 303 mm (a) and annual average precipitation lower than 303 mm (b). The comparison between the calibrated vegetation area extracted from Landsat TM and that extracted from SPOT 5 were checked for regions where annual average precipitation higher than 303 mm (c) and annual average precipitation lower than 303 mm (d). This figure was adapted from the figures of Zhu et al. (2016).

Appendix 4. Assessment of the outcome of the vegetation quality type identified in this study with plot data from National Forest Resource Inventory. Consistency was the ratio of the number of the plots that were identified as healthy vegetation and no degradation, or identified as unhealthy vegetation and severe/potential degradation at the same time to the number of all plots. The assessment was for the outcome of vegetation quality type identification during 2010–2017, for the healthy status data from National Forest Resource Inventory are only available for 2017.

		Vegetation Quality Type
Vegetation Type	Healthy Status	Severe degradation	Potential Degradation A	Potential degradation B	No degradation	Consistency
Forest	Unhealthy Vegetation	51	1	0	0	72.20%
	Healthy Vegetation	6	28	4	255
Shrub-land	Unhealthy Vegetation	8	8	0	111	80.87%
	Healthy Vegetation	74	96	19	1252
Grassland	Unhealthy Vegetation	1	0	0	0	77.14%
	Healthy Vegetation	0	0	7	27
Sum of all vegetation	Unhealthy Vegetation	60	9	0	111	79.62%
	Healthy Vegetation	107	124	30	1507

Appendix 5.

Area of vegetation types of study area in each time point (million ha).

	Vegetation Type	1978	1990	2000	2010	2017
Forest		0.308	0.452	0.508	0.519	0.658
	Coniferous Forest	0.079	0.159	0.180	0.182	0.218
	Broadleaved Forest	0.210	0.252	0.291	0.300	0.410
	Mixed Conifer-broadleaved Forest	0.019	0.041	0.037	0.038	0.030
Shrub-land		0.762	0.937	1.138	2.853	3.344
Grassland		29.867	30.442	29.096	27.720	27.173
Sum of all Vegetation		30.937	31.830	30.742	31.092	31.175

Appendix 6.

TC, NC, and SC of different vegetation types for period of 1978-1990, 1990-2000, 2000-2010, and 2010-2017 (10⁴ ha). TC = total change area; NC = net change area; SC = swap change area.

Period	Vegetation Type	Increase Area	Decrease Area	Total Change Area	Net Change Area	Swap Change Area
1978–1990	Coniferous Forest	8.6	0.6	9.3	8.0	1.3
	Broadleaved Forest	6.8	2.6	9.3	4.2	5.1
	Mixed Forest	2.4	0.2	2.5	2.2	0.3
	Shrub-land	25.1	7.6	32.7	17.5	15.3
1990–2000	Coniferous Forest	3.7	1.6	5.2	2.1	3.1
	Broadleaved Forest	5.0	1.1	6.2	3.9	2.3
	Mixed Forest	0.4	0.8	1.2	−0.4	0.8
	Shrub-land	24.1	4.0	28.1	20.2	7.9
2000–2010	Coniferous Forest	1.6	1.4	3.0	0.2	2.8
	Broadleaved Forest	4.4	3.5	8.0	0.9	7.1
	Mixed Forest	0.3	0.3	0.6	0.0	0.5
	Shrub-land	185.1	13.6	198.7	171.5	27.2
2010–2017	Coniferous Forest	9.6	18.2	27.8	−8.5	19.3
	Broadleaved Forest	28.6	30.0	58.5	−1.4	57.1
	Mixed Forest	1.7	3.8	5.4	−2.1	3.4
	Shrub-land	157.8	285.3	443.1	−127.5	315.5
Summation of 1978–2017	Coniferous Forest	23.5	21.8	45.3	1.7	26.5
	Broadleaved Forest	44.8	37.2	82.0	7.6	71.6
	Mixed Forest	4.8	5.0	9.7	−0.2	5.1
	Shrub-land	392.1	310.5	702.6	81.6	366.0

Appendix 7.

Spatial distribution of each category of vegetation quality types for all vegetation types in four periods during 1978 to 2017.

Appendix 8.

Area of desertification by intensity level during 1978-2017 (million ha).

Desertification Intensity Level	1978	1990	2000	2010	2017
Slight Desertification	6.898	6.321	7.417	8.991	9.050
Moderate Desertification	5.437	6.068	6.017	6.570	6.784
Severe Desertification	4.718	5.999	6.257	6.327	6.351
Extremely Severe Desertification	7.845	8.710	9.480	6.740	6.548
Sum of Desertification	24.898	27.098	29.171	28.628	28.732

Appendix 9. Title and scene ID of the remote sensing images used. For the scene ID, for Landsat images is the “path number and row number”, for SPOT 5 images is the number of “K and J”, and for GF-1 is the “Scene Number”.

Satellite	Scene ID
Landsat TM	118029, 119027, 119028, 119029, 119030, 120027, 120028, 120029, 120030, 121028, 121029, 121030, 121031, 121032, 122026, 122028, 122029, 122030, 122031, 122032, 123025, 123026, 123027, 123028, 123029, 123030, 123031, 123033, 124025, 124026, 124027, 124028, 124029, 124030, 124031, 125025, 125026, 125027, 125029, 125030, 125031,125032, 126029, 126030, 126031, 126032, 126033, 127029, 127030, 127031, 127032, 127033, 128030, 128031, 128032, 128033, 128034, 129031, 129032, 129033, 129034, 130031, 130032, 130033, 130034, 131031, 131032, 131033, 131034, 132031, 132032, 132033, 132034, 132035, 133031, 133032, 133033, 133034 133,035, 134031, 134032, 134033, 134034, 134035, 135032, 135033, 135034, 135035, 136031, 136032, 136033, 136034, 136035, 137031, 137032, 137033, 137034, 137035, 138030, 138031, 138032, 138033, 138034, 139029, 139030, 139032, 139033, 139034, 140029, 140030, 140032, 140033, 140034, 141028, 141029, 141030, 141032, 141033, 141034, 142028, 142029, 142031, 142032, 142033, 142034, 143027, 143029, 143031, 143033, 143034, 144027, 144028, 144029, 144030, 144031, 144032, 144034, 145028, 145029, 145031, 145032, 145033, 145034, 145035, 146029, 146031, 146032, 146033 146,034, 147029, 147031, 147032, 147034, 148032, 148033, 149032, 149033
Landsat MSS	128028, 128029, 128030, 129027, 129028, 129029, 129030, 130028, 130029, 130030, 130031, 131026, 131028, 131029, 131030, 131031, 131032, 132028, 132029, 132030, 132031, 132032, 133026, 133028, 133029, 133030, 133031, 133033, 134026, 134028, 134029, 134030, 134031, 135026, 135029, 135030, 135031, 135032, 136029, 136030, 136031, 137029, 137030, 137031, 137032, 137033, 138030, 138031, 138033, 138034, 139031, 139033,139034, 140031, 140032, 140033, 140034, 141031, 141032, 141033, 141034, 142031, 142032, 142033, 142035, 143031, 143032, 143033, 143034, 143035, 144031, 144032, 144033, 144034, 144035, 145031, 145033, 145034, 145035, 146031, 146032, 146033, 146034, 146035, 147031, 147032, 147033, 147034, 147035, 148030, 148031, 148032, 148033, 148034, 149030, 149032, 149033, 149034, 150030, 150032, 150033, 150034, 151029, 151030, 151032, 151033, 151034, 152029, 152030, 152033, 152034, 153028, 153029, 153031, 153032, 153033, 153034, 154027, 154028, 154029, 154031, 154032, 154034, 155027, 155028, 155029, 155030, 155031, 155032, 155034, 156027, 156028, 156029, 156031, 156032, 156033, 156034, 157027, 157028, 157029, 157031, 157032, 157033, 157034, 158028, 158029, 158031, 158032, 158034, 159031, 159032, 159033, 159035, 160032, 160033, 161032, 161033
SPOT 5	200 267, 200 268, 200 269, 200 270, 203 265, 203 266, 203 267, 203 268, 203 269, 212 252, 212 253, 212 254, 213 252, 213 253, 213 254, 214 254, 214 255, 214 256, 216 252, 216 253, 216 254, 219 254, 219 255, 219 256, 219 257, 219 258, 219 259, 219 260, 219 261, 219 262, 219 263, 219 264, 219 265, 220 260, 220 261, 220 262, 220 263, 220 263, 220 264, 220 265, 221 255, 221 256, 222 265, 222 266, 222 269, 222 270, 222 272, 222 273, 222 274, 223 257, 223 258, 223 259, 223 262, 223 263, 223 264, 223 265, 223 266, 224 259, 224 260, 224 261, 224 262, 224 263, 224 264, 227 265, 227 266, 227 267, 227 271, 227 272, 227 273, 227 274, 227 275, 228 262, 228 263, 228 264, 228 265, 228 265, 228 266, 228 266, 228 267, 228 267, 228 268, 228 272, 228 273, 228 274, 229 264, 229 264, 229 265, 229 266, 229 267, 230 265, 230 266, 230 267, 230 268, 231 253, 231 262, 231 263, 231 264, 231 265, 231 266, 231 267, 231 268, 231 269, 232 262, 232 263, 232 264, 232 272, 232 273, 232 275, 233 263, 233 264, 233 265, 233 266, 233 267, 233 268, 233 276, 233 277, 234 263, 234 264, 234 265, 234 265, 234 269, 234 270, 234 271, 237 276, 237 277, 237 278, 238 264, 238 265, 238 266, 238 267, 238 270, 238 271, 239 273, 239 274, 239 275, 240 275, 240 276, 240 277, 240 278, 242 264, 242 265, 242 266, 242 276, 242 277, 242 278, 242 279, 244 270, 244 271, 244 272, 244 273, 244 274, 244 275, 244 276, 244 277, 244 278, 245 270, 245 271, 245 272, 245 273, 246 264, 246 265, 250 265, 250 267, 250 268, 250 269, 250 270, 251 276, 251 277, 251 278, 251 279, 251 280, 251 281, 252 277, 252 280, 252 281, 253 270, 253 272, 253 273, 253 274, 253 275, 253 277, 254 266, 254 267, 254 268, 254 269, 254 271, 254 273, 255 274, 255 275, 255 276, 255 277, 255 278, 255 279, 255 280, 255 281, 256 274, 256 275, 256 276, 256 277, 257 267, 258 265, 258 266, 258 271, 258 272, 258 273, 258 274, 258 275, 258 276, 260 268, 260 269, 260 270, 260 271, 261 269, 261 270, 261 271, 263 265, 263 266, 263 267, 263 268, 263 270, 263 271, 263 272, 263 273, 263 274, 265 269, 265 270, 265 271, 267 262, 267 263, 267 264, 267 265, 267 266, 267 267, 267 270, 267 271, 267 272, 271 266, 271 267, 271 268, 271 269, 272 258, 272 259, 272 260, 272 261, 272 262, 272 270, 272 271, 274 268, 274 269, 275 268, 278 258, 278 261, 278 262, 279 265, 279 266, 281 258, 281 260, 281 261, 281 262, 281 263, 282 256, 283 255, 283 258, 283 259, 283 260, 283 261, 283 262, 283 263, 284 254, 284 255, 284 256, 284 257, 284 271, 288 255, 288 256, 288 257, 288 258, 293 262, 293 263
GF-1	3795831, 3941683, 3941685, 3941686, 3941699, 3941963, 3941964, 3941965, 3941966, 3941976, 3941977, 3941978, 3941993, 3945318, 3945323, 3945422, 3945438, 3945439, 3947916, 3948032, 3948035, 3949763, 3949764, 3949895, 3949896, 3949897, 3949898, 3952277, 3952289, 3952676, 3952677, 3955954, 3955955, 3955956, 3956214, 3956324, 3956325, 3956326, 3956331, 3956664, 3956665, 3957010, 3957011, 3957012, 3957015, 3957016, 3957017, 3957018, 3957019, 3957021, 3957022, 3957023, 3957024, 3957026, 3957130, 3957131, 3957135, 3957136, 3957137, 3957138, 3957142, 3957144, 3957145, 3959215, 3959216, 3959217, 3959219, 3959220, 3959221, 3959226, 3959345, 3959346, 3959347, 3959348, 3959351, 3959356, 3959358, 3959822, 3966352, 3966353, 3966621, 3966969, 3966970, 3966982, 3967100, 3967101, 3967113, 3967114, 3970785, 3970786, 3970792, 3970793, 3970794, 3970795, 3970808, 3971013, 3971014, 3971019, 3971020, 3971021, 3971022, 3971023, 3971035, 3975227, 3975228, 3975229, 3975354, 3975355, 3975359, 3975832, 3975833, 3975936, 3975937, 3979086, 3979097, 3979098, 3979101, 3979102, 3979215, 3979218, 3979219, 3979220, 3981303, 3981304, 3981306, 3981307, 3981308, 3981309, 3981311, 3981312, 3981478, 3981483, 3981484, 3981485, 3981486, 3981487, 3981488, 3981489, 3981490, 3983073, 3983074, 3983080, 3983081, 3983316, 3983324, 3983325, 3983326, 3983544, 3983545, 3983546, 3983553, 3983554, 3986491, 3986492, 3986495, 3986578, 3986583, 3986595, 3986597, 3986598, 3986599, 3989906, 3989907, 3989908, 3989909, 3989910, 3989913, 3989914, 3990096, 3990098, 3990105, 3990107, 3990202, 3990208, 3990211, 3993149, 3993150, 3993151, 3993613, 3993618, 3993630, 3993631, 3993632, 3993634, 3993638, 3993639, 3993646, 3993732, 3993747, 3993748, 3993749, 3993750, 3993751, 3993752, 3993753, 3993754, 3993755, 3993756, 3997094, 3997095, 3997098, 3997319, 3997320, 3997321, 3997324, 3997325, 3997338, 3997552, 3997555, 3997556, 3997558, 3997559, 3997560, 3997570, 4001116, 4001117, 4001118, 4002599, 4002600, 4002601, 4002602, 4002856, 4002857, 4002858, 4002859, 4005740, 4005741, 4005742, 4005743, 4005745, 4005746, 4005747, 4005748, 4005756, 4005757, 4005758, 4005759, 4005762, 4005763, 4005766, 4005768, 4005774, 4005775, 4006025, 4006026, 4006027, 4006029, 4006031, 4006032, 4006033, 4006041, 4006043, 4006045, 4006049, 4006050, 4006052, 4006053, 4006054, 4006167, 4006265, 4006266, 4006273, 4009213, 4009217, 4009218, 4009327, 4009328, 4009330, 4009331, 4009597, 4009598, 4009599, 4009689, 4013730, 4013732, 4014152, 4014153, 4014164, 4014165, 4014166, 4016867, 4016868, 4016869, 4016872, 4016874, 4017310, 4017312, 4017314, 4017316, 4017318, 4017420, 4017421, 4020727, 4020728, 4020737, 4020738, 4020739, 4020742, 4020743, 4020744, 4020755, 4020980, 4020981, 4020982, 4020991, 4020992, 4020993, 4020994,
	4020996, 4020998, 4021008, 4021090, 4021091, 4021099, 4021309, 4021310, 4021318, 4024359, 4024360, 4024361, 4024666, 4024667, 4024669, 4027303, 4027304, 4027310, 4027311, 4027406, 4027408, 4027412, 4027414, 4027415, 4027416, 4033425, 4036282, 4036283, 4036284, 4036285, 4036295, 4036511, 4036512, 4036513, 4036514, 4036515, 4036516, 4036517, 4036518, 4036519, 4036811, 4036812, 4036813, 4036814, 4036931, 4036932, 4036933, 4038862, 4038863, 4039243, 4039244, 4039485, 4039486, 4039487, 4039489, 4039490, 4039603, 4039604, 4039606, 4039607, 4039608, 4039715, 4039716, 4039717, 4039718, 4039851, 4039853, 4044028, 4044029, 4044030, 4044031, 4044033, 4044119, 4044120, 4044121, 4044122, 4044123, 4044124, 4044125, 4044126, 4044131, 4044422, 4044423, 4044430, 4044431, 4044443, 4044445, 4044775, 4044776, 4044777, 4044781, 4044784, 4044796, 4044797, 4044798, 4044810, 4044812, 4046993, 4046994, 4046996, 4046997, 4046998, 4046999, 4047000, 4047002, 4047003, 4047004, 4047005, 4047006, 4047008, 4047378, 4047379, 4047380, 4047390, 4047473, 4047475, 4051180, 4051182, 4051183, 4051186, 4051187, 4051193, 4051194, 4051196, 4051198, 4051199, 4051202, 4051205, 4051206, 4051462, 4051463, 4051464, 4051465, 4051466, 4051469, 4051474, 4051475, 4051478, 4051479, 4051480, 4051481, 4051482, 4051483, 4051484, 4051485, 4051486, 4051487, 4051489, 4052148, 4052259, 4052260, 4055279, 4055394, 4055396, 4055397, 4055398, 4055851, 4055852, 4055853, 4055854, 4055855, 4055856, 4055857, 4055869, 4055961, 4055968, 4055969, 4055972, 4055973, 4055975, 4055976, 4055977, 4055987, 4059548, 4059549, 4059560, 4059638, 4059639, 4059640, 4059641, 4059642, 4059651, 4059652, 4059854, 4059858, 4059861, 4059863, 4060078, 4060079, 4060080, 4060082, 4060083, 4060084, 4060085, 4060086, 4065553, 4065640, 4065669, 4065670, 4065671, 4065673, 4065675, 4066126, 4066127, 4066128, 4066136, 4066140, 4066141, 4066142, 4066143, 4066144, 4066146, 4066148, 4066149, 4066308, 4066309, 4066316, 4066317, 4066318, 4066323, 4066324, 4066326, 4066327, 4066329, 4066330, 4066332, 4066759, 4066760, 4066769, 4066770, 4066947, 4066949, 4066956, 4066957, 4071387, 4071389, 4071390, 4071392, 4071394, 4071396, 4071397, 4071398, 4071399, 4071401, 4071402, 4071403, 4071404, 4071626, 4071628, 4071630, 4071631, 4071632, 4071633, 4071639, 4071640, 4071641, 4071719, 4071721, 4071733, 4071734, 4071735, 4071837, 4071838, 4071839, 4071851, 4071853, 4071863, 4071875, 4075598, 4075601, 4075602, 4075603, 4075604, 4075605, 4075607, 4075608, 4075610, 4075702, 4075703, 4075704, 4075705, 4075706, 4075707, 4075709, 4075713, 4075716, 4075717, 4075718, 4076369, 4076372, 4078936, 4078937, 4078938, 4078940, 4078941, 4078942, 4078944, 4078945, 4078946, 4078948, 4078949, 4079056, 4079057, 4079058, 4079059, 4079060, 4079061, 4079062, 4079063, 4079064, 4079065, 4079067, 4079375, 4079376, 4079434, 4079438, 4082994, 4083100, 4084339, 4084341, 4084343, 4084344, 4084345, 4084357, 4084358, 4084359, 4084361, 4084362, 4084363, 4084364, 4084365, 4084366, 4084459, 4084460, 4084461, 4084462, 4084468, 4084469, 4084470, 4084471, 4084472, 4084473, 4084475, 4084478, 4084479, 4084481, 4084485, 4086142, 4086331, 4086332, 4086341, 4086343, 4086344, 4088164, 4088165, 4088293, 4088303, 4088304, 4088828, 4091621, 4091622, 4091780, 4091784, 4091785, 4091786, 4092080, 4092082, 4092095, 4092096, 4092097, 4092110, 4092322, 4092323, 4092328, 4092330, 4092331, 4092344, 4092345, 4092359, 4094588, 4094589, 4094590, 4094593, 4094594, 4094596, 4094652, 4094653, 4094654, 4094656, 4094657, 4094659, 4094660, 4094661, 4096906, 4096907, 4096909, 4096910, 4096911, 4096912, 4096914, 4096915, 4096916, 4096917, 4096919, 4097038, 4097039, 4097040, 4097041, 4097042, 4097044, 4097045, 4097046, 4097047, 4097050, 4097051, 4097052, 4097289, 4097291, 4097298, 4097299, 4097300, 4097376, 4097377, 4097378, 4097385, 4097386, 4100656, 4100657, 4100658, 4100660, 4100661, 4100662, 4100671, 4100673, 4100851, 4100863, 4100864, 4101118, 4101139, 4101141, 4101142, 4104059, 4104060, 4104061, 4104062, 4104064, 4104076, 4104077, 4105247, 4105248, 4105252, 4105253, 4105351, 4105353, 4105354, 4105355, 4105518, 4105522, 4105523, 4105525, 4105526, 4105536, 4105537, 4105777, 4105782, 4105783, 4105784, 4105786, 4105787, 4105789, 4105790, 4105791, 4105792, 4105802, 4111011, 4111016, 4111030, 4112100, 4112101, 4112102, 4112103, 4112104, 4112104, 4112106, 4112106, 4114688, 4114689, 4114693, 4115143, 4115144, 4115146, 4115147, 4115148, 4115149, 4115150, 4115151, 4115157, 4115398, 4115399, 4115400, 4115402, 4115403, 4115404, 4115405, 4115412, 4116560, 4116561, 4116563, 4116564, 4116566, 4118372, 4118389, 4118390, 4118461, 4118462, 4118465, 4118466, 4118477, 4118478, 4118479, 4118496, 4119094, 4119107, 4119352, 4119364, 4119365, 4119374, 4121192, 4121194, 4121195, 4121196, 4121198, 4121199, 4121200, 4121201, 4121202, 4121203, 4121204, 4121206, 4121208, 4121209, 4121441, 4121442, 4121444, 4121445, 4121446, 4121450, 4121451, 4121453, 4121454, 4121455, 4121456, 4121460, 4121461, 4121462, 4121463, 4123340, 4123341, 4123342, 4123343, 4123344, 4123345, 4123346, 4123348, 4123549, 4123551, 4123553, 4123554, 4123555, 4123556, 4123557, 4124145, 4126033, 4126037, 4126039, 4126040, 4126049, 4126050, 4126052, 4126053, 4129609, 4129610, 4129611, 4129613, 4129614, 4129616, 4129617, 4129619, 4129862, 4129863, 4129864, 4129865, 4129866, 4129867, 4129868, 4129869, 4129870, 4129871, 4129872, 4129873, 4129874, 4130144, 4130145, 4130248, 4130256, 4132895, 4132896, 4132899, 4132900, 4132901, 4132902, 4132903, 4132910, 4132911, 4132912, 4132914, 4132915, 4134988, 4134990, 4134991, 4134994, 4135004, 4135005, 4135006, 4135008, 4135009, 4135011, 4135012, 4135014, 4135015, 4135017, 4135022, 4135109, 4135110, 4135123, 4135124, 4135125, 4135126, 4135127, 4135128, 4135129, 4135130, 4135131, 4135132, 4135141, 4135143, 4135277, 4135278, 4135279, 4135298, 4135444, 4135446, 4135447, 4135466, 4139141, 4139143, 4139144, 4139146, 4139148, 4139149, 4139151, 4139152, 4139390, 4139392, 4139395, 4139396, 4139397, 4139398, 4139399, 4139400, 4139403, 4139404, 4139661, 4139662, 4139663, 4139664, 4139674, 4139675, 4139677, 4139678, 4139688, 4139781, 4139794, 4139795, 4139796, 4139797, 4139805, 4141446, 4141447, 4141448, 4143355, 4143356, 4143357, 4143358, 4143359, 4143363, 4143364, 4143365, 4143378, 4143461, 4143462, 4143465, 4143466, 4143473, 4143481, 4143482, 4143483, 4143484, 4143485, 4143843, 4143844, 4143846, 4143847, 4143848, 4143849, 4146742, 4146743, 4146744, 4146745, 4146746, 4146747, 4146748, 4146750, 4146751, 4146752, 4146834, 4146834, 4146835, 4146836, 4146837, 4146838, 4146840, 4146841, 4146843, 4146844, 4146845, 4146847, 4151499, 4151500, 4151502, 4151505, 4151506, 4151507, 4151508, 4151509, 4151510, 4151512, 4151517, 4151585, 4151586, 4151587, 4151588, 4151589, 4151590, 4151591, 4151601, 4151602, 4151603, 4151604, 4151605, 4151607, 4151608, 4151611, 4151612, 4151618, 4151619, 4151892, 4151894, 4151915, 4152137, 4152138, 4152140, 4152161, 4153344, 4153345, 4153347, 4153348, 4153357, 4153358, 4153584, 4153585, 4153587, 4153596, 4154079, 4154182, 4154197