LISA is a local statistical method for spatial variables that reveals the spatial clustering characteristics of observations within spatially adjacent regions [33,60]. A positive LISA value indicates that the SDI is similar to the adjacent value and reveals a spatial pattern of high–high clustering (H–H, high values are near other high values) or low–low clustering (L–L, low values are near other low values). A negative LISA value indicates that the SDI is a spatial outlier and can include a high–low outlier value (H–L, a high value is near a low value) and a low–high outlier value (L–H, a low value is near a high value).

#### **3. Results**

#### *3.1. Integration of the Remote Sensing Indexes Based on PCA*

The SDI of the ADD was calculated by the PC1 of the four indicators. The PC1 results shown in Table 1 indicated that, during the studied years, the percent eigenvalues of PC1 were higher than 78%, revealing that PC1 integrated most of the information of the four indictors. Therefore, PC1 was chosen to construct the SDI in this study. The loading values of the four variables in PC1 were divided into two types according to their signs. The albedo and SI comprised one type with loading values that were positive, and the NDVI and LSM comprised the second type with loading values that were negative. The opposite

signs of the two variables indicated that the corresponding contributions to the SDI value were opposing.

**Table 1.** Loadings of the four selected variables on the first principal component (PC1) and associated contributions in different study years.

SI: salinization index (SI); NDVI: normalised difference vegetation index; LSM: land surface soil moisture index.

The descriptive statistics of the PC1 indicated that the average PC1 value increased from 0.42 in 2010 to 0.43 in 2000, and the medium value of PC1 also increased from 0.38 to 0.40 (Table 2). From 2010 to 2019, on the contrary, the average PC1 value decreased from 0.41 to 0.30, and the medium value of PC1 also decreased from 0.38 to 0.28. The positive skewness in 1990, 2000, 2010 and 2019 indicated that the tail on the right side of the probability density function was longer or fatter than the left side. Overall, the statistical results of PC1 showed that land degradation accelerated from 1990 to 2000 and weakened from 2010 to 2019.

Furthermore, the correlation coefficient between each indicator and the SDI and that among the indicators are shown in Figure 3 (at the 0.01 level of significance). For four years, the SDI exhibited a high correlation with each single indicator. In general, the SDI exhibited a positive correlation with the SI and albedo and a negative correlation with the NDVI and LSM (Figure 3). The SI exhibited the highest correlation with the SDI, and the positive correlations were 0.971, 0.983, 0.973 and 0.972 in 1990, 2000, 2010 and 2019, respectively. The albedo exhibited the highest negative correlation with the SDI in 2000, and the correlation was 0.915. The NDVI exhibited the highest correlation with the SDI in 2019 (−0.885). The correlation coefficients between the LSM and SDI were greater than −0.72, and the highest correlation with the SDI was observed in 2000 (−0.911).

**Figure 3.** Correlations between pairs of the four selected indicators and their correlations with the SDI in different study years. SI: salinization index; NDVI: normalised difference vegetation index; LSM: land surface soil moisture index; SDI: salinized land degradation index. The blue and red colours represent negative and positive correlations, respectively (the darker the colour, the stronger the correlation).

#### *3.2. Spatiotemporal Changes in the Land Degradation*

To analyse the spatiotemporal characteristics of the land degradation during the different periods in the ADD, the SDI values were normalised by Equation (6) (range of 0 to 1). As the SDI approximates a normal distribution, we divided it into five categories by equal intervals to indicate the different land degradation levels [36,58]—namely, no degradation (0–0.2), slight degradation (0.20–0.4), moderate degradation (0.4–0.6), strong degradation (0.6–0.8) and extreme degradation (0.8–1). In summary, the land degradation level distribution was not uniform in space and varied over space and time. As shown in Figure 4, in terms of the land degradation level distribution, the extreme and strong degradation areas were clustered in the west and north of the ADD during the studied years. Areas with moderate degradation corresponded to a sporadic distribution in the middle and south of the ADD. Most areas in the middle of the ADD exhibited a slight degradation. The spatial distribution of the not-degraded areas in the study area showed a difference in the 4 years: in 1990 and 2000, the not-degraded areas were mainly distributed along the Amu Darya River, and a small portion appeared in the northwest and northeast corners of the ADD; in contrast, the not-degraded areas were mainly distributed in the middle of the study area in 2010 and 2019 and formed a "V" shape.

**Figure 4.** Spatial distribution of the land degradation levels in the ADD in each study year. Extreme: extreme degradation; Strong: strong degradation; Moderate: moderate degradation; Slight: slight degradation.

Figure 5 shows the percentage of the study area occupied by the five land degradation levels in 1990, 2000, 2010 and 2019. In general, the largest areas in the ADD were slight degradations in the years 1990–2019, which accounted for more than 26% of the total area covered by the Landsat images. From 1990 to 2000, the areas with extreme, strong and no degradation decreased, whereas the areas with slight and moderate degradations increased. From 2000 to 2010, the areas with extreme, strong and no degradation expanded; among which, the expansion of the no degradation regions was significant. In contrast, the areas

with slight and moderate degradations exhibited a decreasing trend; in particular, for the areas with a slight degradation, the dynamic degree was 7%. From 2010 to 2019, a small increase from 13.30% to 15.56% of the total area was observed in the area with moderate degradation, and the areas of the other four levels did not change considerably.

**Figure 5.** Proportion of land degradation levels in different study years. Extreme: extreme degradation; Strong: strong degradation; Moderate: moderate degradation; Slight: slight degradation.

Using a spatial analysis, the land degradation spatial distribution changes from 1990 to 2000, 2000 to 2010, 2010 to 2019 and 1990 to 2019 were mapped (Figure 6). We defined the figure elements as follows: development of land degradation across levels 1 or 2 corresponded to "Developed" (e.g., a change from "No degradation" to "Slight degradation" or "Moderate degradation"), a development across levels 3 or 4 corresponded to "Seriously developed" (e.g., a change from "No degradation" to "Strong degradation" or "Extreme degradation"), an improvement of the land degradation across levels 1 or 2 corresponded to "Improvement" (e.g., a change from "Extreme degradation" to "Strong degradation" or "Moderate degradation"), an improvement across levels 3 or 4 corresponded to "Significant improvement" (e.g., a change from "Extreme degradation" to "Slight degradation" or "No degradation"), and no change during the study periods corresponded to "Stable". The areas of land degradation dynamics for the four periods are shown in Table 3.

**Table 3.** Area and proportion of land degradation developed or an improvement of the ADD during different time periods.

**Figure 6.** Spatial distribution of land degradation developed or improvement in the ADD in different time periods.

Overall, the stable areas accounted for a large proportion (more than 42%) during the four study periods (Table 3), and these areas were mainly located west and southeast of the ADD. From 1990 to 2000, the improvement areas covered 3653.78 km2 (30.6%), and these areas were mainly clustered north of the ADD (Figure 6). Additionally, the developed areas covered 2687.53 km2 (22.5%), and these areas mainly occupied the south and middle regions of the ADD. A smaller proportion (1.4%) of significant improvement areas was observed in the northern part of the ADD. In comparison, the proportion of the seriously developed areas was smaller (1.1%), and these areas were mainly located in the northwest corner of the ADD. From 2000 to 2010, the developed areas covered 3728.76 km2 (31.3%), and the areas were mainly concentrated in the northern part of the study area. The improvement areas covered 2886.02 km2 (24.2%) from 2000 to 2010 and were mainly located in the eastern and central parts of the ADD. The areas of seriously developed and significant improvement exhibited smaller proportions—0.92% and 0.88%, respectively—and the seriously developed areas were mainly observed downstream of the ADD. Compared with those in 2000–2010, during 2010–2019, the areas with developed, and the improvements exhibited decreasing trends and occupied 20.5% and 19.7% of the total area, respectively. Moreover, these regions were mainly clustered in the downstream of the ADD. The seriously developed and significant improvement areas occupied less than 0.4% of the total area. From 1990 to 2019, the developed areas covered 2901.65 km2 (24.3%) and were mainly concentrated north of the ADD and in the downstream region of the Amu Darya River. The improvement regions, which occupied 27.8% of the total area, were mainly observed in the west and east of the ADD. The seriously developed regions, with an area of approximately 217.75 km2, were mainly concentrated in the north of the ADD, and the areas with significant improvement were smaller, accounting for only 0.9% of the total area.

#### *3.3. Spatial Autocorrelation Analysis of the SDI*

To further clarify the spatial and temporal variabilities in the land degradation, the spatial autocorrelation of the SDI was examined.

The mapping of the global spatial autocorrelation of the SDI is shown in Figure 7. Most of the SDI values are distributed in the first and third quadrants, with H–H and L–L clustered in the first and third quadrants, respectively, indicating a strong positive spatial correlation between the spatial units in these two quadrants. The Moran's I values in 1990, 2000, 2010 and 2019 were 0.888, 0.856, 0.891 and 0.851, respectively, which were high values greater than zero. The results showed that the SDI of the ADD exhibited significant spatial clustering, which indicated a strong positive spatial correlation. The Moran's I values of the SDI decreased, increased and later decreased again in 1990, 2000, 2010 and 2019, exhibiting an overall decreasing trend during the study period.

**Figure 7.** Moran scatter plot of the SDI in the ADD for each study year.

*3.4. Spatial and Temporal Changes in Land Use and Salinization*

Figure 8 displays the land use maps for 1990, 2000, 2010 and 2019. These land use maps indicated that croplands were the dominant land use type in the ADD. Grassland and forest were distributed in the north of the ADD and bare soil mainly in the edge and north of the ADD. Figure 8b displays the spatial variations in the land use types. The combinations with no land use type transformations, smaller conversion areas and built-up land transformations were merged into "Stable and others". From 1990 to 2010, the conversion of land use categories was mainly between grasslands and croplands. The conversion from grassland to cropland (GL to CL) was prominent in the northern part of the ADD, where the area of grassland decreased by 158.74 km2 and 338.93 km2 from 1990 to 2000 and 2000 to 2010, respectively (Table 4). However, the conversion of land use categories from 2010 to 2019 was mainly cropland to grassland (CL to GL), which was distributed in the northern and western parts of the ADD. The area of cropland decreased by 602.26 km2, and the area of grassland increased by 492.37 km2 during this period. In addition, there was partial degradation of grassland to bare soil (GL to BS) in the northern

part of the ADD from 2010 to 2019. The conversion of land use categories throughout the study period was mainly from grassland to cropland and cropland to grassland.

**Figure 8.** Land use change maps from 1990 to 2019. CL: cropland; GL: grassland; FR: forest; BS: bare soil. (**a**) land use maps for 1990, 2000, 2010 and 2019; (**b**) changes between different periods.

**Table 4.** Areas and percentages of different land use types from 1990 to 2019 in the ADD.

Figure 9 reveals the spatial distribution (Figure 9a) and variation (Figure 9b) of the SI over the study time period. The values of the SI were larger at the edges and north of the ADD and smaller in the centre. The SI in the northern part of the ADD decreased from 1990 to 2000, while it increased significantly from 2000 to 2010. From 2010 to 2019, the SI increased slightly in the central part of the ADD and decreased in the northeast. Throughout the study period (1990–2019), the SI increased significantly in the northern part of the ADD and decreased at the edges. To reveal the relationship between land use change and salinization, the mean SI values for each land use type (built-up land and water were excluded) were calculated in Figure 10. The mean value of the SI for the land use type was the largest in 2000 compared to the other three years, indicating higher soil salinization. The mean value of the SI for the land use type was the largest in 2000 compared to the other three years, indicating higher soil salinization, with a decreasing trend in the SI from 2000 to 2019. Within each year, bare soil had the largest mean value of the SI, followed by grassland and forest.

The dynamic characteristics of salinization during land use change were illustrated in Figure 11. The ΔSI represents the difference in the SI over the study period (the following year minus the previous year). We counted a percentage of areas with ΔSI < 0 and ΔSI < 0 over the course of each land use type transfer. The SI increased during all land use type changes from 1990 to 2000 (more than 50% of the area with ΔSI > 0), which was related to the maximum SI value in 2000 mentioned earlier. The 2000–2010 SI values decreased for BS-CL and GL-CL (ΔSI < 0 over 80% of the area). The area with ΔSI < 0 dominated the land use change process from 2010 to 2019. During the conversion of bare soil to cropland (BS-CL), forest (BS-FR) and grassland (BS-GL) from 1990 to 2019, the area with ΔSI < 0 exceeded 50%, indicating a decrease in the SI.

**Figure 9.** Spatial distribution of soil salinization from 1990 to 2019 (**a**) and spatial distribution of the salinity index (SI) changes between different periods (**b**).

**Figure 10.** Mean values of the salinity index (SI) for land use types over the study period.

**Figure 11.** The percentage of ΔSI area for the land use categories that changed during the study period. ΔSI: the salinity index (SI) of the next year minus the SI of the previous year. BS–CL: bare soil to cropland; BS–FR: bare soil to forest; CL–GL: cropland to grassland; FR–BS: forest to bare soil; GL–BS: grassland to bare soil; GL–CL: grassland to cropland; BS–GL: bare soil to grassland.

#### **4. Discussion**

#### *4.1. Effectiveness of the Proposed SDI*

In this study, we established a new SDI by integrating the SI, albedo, NDVI and LSM indices based on the PCA method. The SDI was used to explore the land degradation characteristics for a typical salinized area, i.e., the ADD. We found that the regions with extreme land degradation were mainly distributed downstream and at the periphery of the ADD (Figure 4). The SDI-based results supported previous findings that the ecological risks and vulnerabilities are higher in the downstream and peripheral regions of the ADD [10,46], thereby demonstrating that the SDI can reflect the land degradation conditions of the ADD. To further evaluate the effectiveness of the proposed SDI, we evaluated its effectiveness by field survey data (Figure S1 and Table S1). The relationship between the SDI and soil salt content shown in Figure 12 indicates that the soil salt content was significantly positively correlated with the SDI (R2 = 0.89, *p* < 0.001). The relationships between the four indices derived from Landsat imagery in 2019 and field-measured soil salinity are also presented in Figure S2. There was a positive correlation between the SI, albedo and measured soil salinity (R2 = 0.41 and 0.43, respectively) and a negative correlation between the NDVI, LSM and measured soil salinity (R2 = 0.29 and 0.19, respectively). In contrast to the single indices, the SDI showed a stronger agreement with the measured soil salinity. In general, although the land degradation is influenced by multiple aspects of the environment, this positive correlation suggests that the SDI can capture the salinity features pertaining to the land degradation, which provides potential evidence for the effectiveness of the index in monitoring the land degradation in salinized areas. In addition, the SI extracted from the remote sensing data exhibits a positive correlation with the SDI (Figure 3). The reliability of the SDI is also reflected in the other three indicators. High vegetation cover and sufficient soil moisture reduce the risk of land degradation, and this finding is supported by the negative correlation between the SDI and NDVI and LSM in our study (Figure 3). An increase in the albedo values leads to a higher SDI, which is related to the exposure

information represented by the albedo and may be attributed to a strong coupling between the soil salinity information and albedo (Table 1 and Figure 3). These characteristics have also been reported previously [33]. These results indicate that the SDI can help reliably and efficiently monitor the land degradation of salinized areas. In addition, the SDI is composed of accessible remote sensing indicators and can thus be extended to other similar ecological environments [35,59,63].

**Figure 12.** Relationship between the SDI and field-measured soil salt content.

#### *4.2. Factors Influencing the Land Degradation in the ADD*

Climate variables mainly affect the land degradation through changes in the precipitation and temperature [10,64]. The annual mean precipitation (AMP) and annual mean temperature (AMT) of Nurkus Weather Station are shown in Figure 13. It can be noted that the AMP decreased and AMT increased in the ADD in the past 40 years. From 1990 to 2019, land degradation developed over an area of more than 3000 km2 (Table 3). Previous studies have raised concerns regarding the withering of grasslands and sparse vegetation caused by warming and dry climates in the ADD, warning that these aspects could accelerate land degradation [10,65].

**Figure 13.** Changes in the annual mean precipitation and temperature measured at Nurkus Station from 1980 to 2016.

Amu Darya River is the main source of water for ADD living and irrigation. Overwatering has disrupted the water system of the ADD, leading to the disappearance of an amount of lakes, followed by local climate change, which has reduced the ecosystem stability of the ADD, particularly in the downstream region [7,10]. The widespread use of diffuse irrigation has caused individuals to compensate for inefficient irrigation by collecting large amounts of water, resulting in a reduction in the ecological water that sustains the ecosystems and the gradual withering of vegetation without sufficient water to support growth [9,66,67], exacerbating the land degradation. The extent of the land degradation followed the same trend as that of the water withdrawal from Amu Darya River. With the decline in the water withdrawal in 1990–2000 (Figure 14), the improvement in ADD land degradation was most pronounced, while, in 2000–2019, as the water withdrawal increased, the land degradation developed in a larger area than the improvement area. A more critical situation is that the reservoir built in the upper reaches of the Amu Darya River intercepted a large amount of the water [66], resulting in a decrease in the supply of ecological water downstream. The ecological effects caused by these factors were confirmed by our research: The land degradation downstream of the ADD was significantly degraded compared with the region in the study period (Figures 4 and 6). In addition, the higher levels of land degradation in the outer delta, farther from Amu Darya River, were likely caused by the lack of water supply to the ecosystem and the difficulties in the land management in the transboundary area [46,68]. This finding indicates that, to alleviate the ADD land degradation and the ecological crisis in the Aral Sea basin, further effort and cooperation is necessary in the rational allocation of water resources.

**Figure 14.** Annual changes in the water withdrawal and salt discharge in the ADD from 1990 to 2015.

In addition, the impacts of land-use changes on the land degradation cannot be ignored [10,69,70]. The consolidation and management of croplands can contribute to the mitigation of land degradation, and our research supports this perspective. We demonstrated that the area with no degradation occupied a larger proportion of cropland during most of the study period, while the area with extreme land degradation was mainly distributed on bare soil (Figure 15). In general, the risk of land degradation was reduced when land with sparse vegetation and bare soil was reclaimed as cropland, as crops contribute to higher ecosystem productivity and stability. Land degradation is more severe in the northern part of the ADD, where part of the cropland has been abandoned and converted into grassland or bare soil (Figure 8). Previous studies have indicated that the ADD is facing an ecological threat posed by the degradation of grasslands and croplands to bare soil [10,46]. Compared to the other land use categories, bare soil has a higher soil salinity

(Figure 10), and the conversion of land use types to bare soil not only reduces the biomass but also increases the risk of soil salinization.

**Figure 15.** Percentage graphs showing the proportion of each land degradation level in different land use types. Extreme: extreme degradation; Strong: strong degradation; Moderate: moderate degradation; Slight: slight degradation.CL: cropland; FR: forest; GL: grassland; BS: bare soil.

However, to reduce the salt content of croplands, a large amount of water is acquired from Amu Darya River to rinse the cropland soil, which further aggravates the water deficit of other ecosystems. In addition, the excess salinity from croplands is discharged by widely distributed channels to the Amu Darya River, as well as to the lakes downstream of the delta, resulting in increased salinity in the river water and a significant accumulation of salinity downstream (Figure 9) [10]. The increase in the discharge was the most significant after 2000 (Figure 14). The extreme land degradation distribution patterns were noted to be clustered downstream of the ADD (Figures 4 and 6). These ecological effects caused by the large accumulation of salt in the downstream region were confirmed by our research. In addition, with the disintegration of the Soviet Union, the socialist economy turned into a market economy, and the gradual influx of the rural population into cities led to certain croplands eventually transforming into unused lands with a low biodiversity [71], thereby accelerating land degradation.

The 15th initiative of the Sustainable Development Goals (SDGs) aims to achieve land degradation neutrality by 2030. To effectively alleviate land degradation of the ADD and promote the further realisation of SDGs, based on the abovementioned factors influencing the SDI, the following corresponding measures and countermeasures are proposed. First, salinization treatment technology should be implemented to alleviate the promotion of land degradation caused by salinity, and reservoirs for storing salt can be built to reduce the ecological pressure caused by the transportation of salt from the alkali drainage canal that discharges into the downstream area of the ADD. Second, a drip irrigation system can be promoted to achieve precision irrigation and enhance the irrigation efficiency to relieve

the pressure of water resources required to maintain the stability of the land degradation. Furthermore, farmers are encouraged to maintain the stability and biodiversity of croplands through agricultural subsidy policies. In addition, ecological conservation projects can be considered to mitigate the impacts of climate change on the land degradation.

#### **5. Conclusions**

We coupled multiple remote sensing indices (SI, NDVI, albedo and LSM) to construct a new SDI by using the PCA method. The proposed index integrated the soil salinity, soil bareness, soil moisture and vegetation coverage and made it possible to identify the characteristics of the regional land degradation, especially in salinized areas. To test the reliability of the SDI, the index was applied to the typical ADD region to monitor the spatial and temporal dynamics of the land degradation.

The results indicated that the NDVI and LSM adversely influenced the land degradation, while the SI and albedo had positive effects. The SI was strongly positive correlated with the SDI, with an average correlation coefficient of 0.97. Regions with extreme and strong land degradation were mostly clustered west and north of the ADD. The temporal and spatial dynamics of the SDI indicated that the land degradation in the ADD developed by approximately 26% (including seriously developed and developed areas) from 1990 to 2019, and the degradation was mainly concentrated in the downstream region of the ADD. The areas exhibiting improvement accounted for approximately 28% of the total area of the ADD and were mainly centred in the eastern and central parts. Among them, the area of land degradation developed from 2000 to 2010 was the largest (approximately 32%), while the improvement area was 25%. The results of spatial autocorrelation analysis showed that the SDI values of Moran's I in 1990, 2000, 2010 and 2019 were 0.89, 0.86, 0.89 and 0.85, respectively, which showed that the SDI was clearly clustered in space rather than randomly distributed.

The drying climate and excessive water withdrawal from the Amu Darya River exacerbated the land degradation in the ADD, especially in 2000–2019; as the water withdrawal increased, the land degradation developed into a larger area than the improvement area. The expansion of unused land increases the risk of land degradation, with higher levels of land degradation on unused land than on other types during the study period. In addition, a large amount of salt discharged from croplands downstream of the ADD results in the downstream being the most degraded area of land.

Clarifying the characteristics of land degradation of salinized areas is conducive to the restoration and promotion of sustainable terrestrial ecosystems. Our study revealed land degradation characteristics at the interannual scale of the ADD based on the SDI, which provided an efficient decision-making basis for regional land management. Nevertheless, some limitations still exist in this research. For example, the seasonal and continuous dynamics of land degradation have not been taken well into consideration due to the limited temporal resolution of the Landsat satellites. Constructing a SDI with the high temporal resolution MODerate resolution Imaging Spectroradiometer (MODIS) has the potential to enable the seasonal and continuous temporal monitoring of land degradation on a large scale, which is further work that deserves to be advanced.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/rs13152851/s1: Figure S1: Soil sampling sites of the ADD. Table S1: Field sampling data of soil salt.

**Author Contributions:** T.Y. and G.J. designed the research. T.Y. processed the data and wrote the manuscript. A.B. and G.Z. revised the manuscript. L.J., Y.Y. and X.H. provided the analysis tools and technical assistance. All authors contributed to the final version of the manuscript by proofreading and offering constructive comments. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDA19030301) and the Open Foundation of State Key Laboratory of Desert and Oasis Ecology, Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences (G2019-02-03).

**Data Availability Statement:** The data is available upon request.

**Acknowledgments:** We thank the journal's editors and reviewers for their kind comments and valuable suggestions to improve the quality of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Soil Salinity Assessment in Irrigated Paddy Fields of the Niger Valley Using a Four-Year Time Series of Sentinel-2 Satellite Images**

**Issaka Moussa 1,2, Christian Walter 1,\*, Didier Michot 1, Issifou Adam Boukary 3, Hervé Nicolas 1, Pascal Pichelin <sup>1</sup> and Yadji Guéro <sup>2</sup>**

Received: 31 August 2020; Accepted: 8 October 2020; Published: 16 October 2020

**Abstract:** Salinization is a major soil degradation threat in irrigated systems worldwide. Irrigated systems in the Niger River basin are also affected by salinity, but its spatial distribution and intensity are not currently known. The aim of this study was to develop a method to detect salt-affected soils in irrigated systems. Two complementary approaches were tested: salinity assessment of bare soils using a salinity index (SI) and monitoring of indirect effects of salinity on rice growth using temporal series of a vegetation index (NDVI). The study area was located south of Niamey (Niger) in two irrigated systems of rice paddy fields that cover 6.5 km2. We used remote-sensing and ground-truth data to relate vegetation behavior and reflectance to soil characteristics. We explored all existing Sentinel-2 images from January 2016 to December 2019 and selected cloud-free images on 157 dates that covered eight successive rice-growing seasons. In the dry season of 2019, we also sampled 44 rice fields, collecting 147 biomass samples and 180 topsoil samples from January to June. For each field and growing season, time-integrated NDVI (TI-NDVI) was estimated, and the SI was calculated for dates on which bare soil conditions (NDVI < 0.21) prevailed. Results showed that since there were few periods of bare soil, SI could not differentiate salinity classes. In contrast, the high temporal resolution of Sentinel-2 images enabled us to describe rice-growing conditions over time. In 2019, TI-NDVI and crop yields were strongly correlated (r = 0.77 with total biomass yield and 0.82 with grain yield), while soil electrical conductivity was negatively correlated with both TI-NDVI (r = −0.38) and crop yield (r = −0.23 with total biomass and r = −0.29 with grain yield). Considering the TI-NDVI data from 2016–2019, principal component analysis followed by ascending hierarchical classification identified a typology of five clusters with different patterns of TI-NDVI during the eight growing seasons. When applied to the entire study area, this classification clearly identified the extreme classes (i.e., areas with high or no salinity). Other classes with low TI-NDVI (i.e., during dry seasons) may be related to areas with moderate or seasonal soil salinity. Finally, the high temporal resolution of Sentinel-2 images enabled us to detect stresses on vegetation that occurred repeatedly over the growing seasons, which may be good indicators of soil constraints due to salinity in the context of the irrigated paddy systems of Niger. Further research will validate the ability of the method developed to detect moderate soil salinity constraints over large areas.

**Keywords:** salinization; irrigated systems; Niger River basin; salinity index; vegetation index; TI-NDVI; Sentinel-2 images; high temporal resolution

#### **1. Introduction**

Irrigated agriculture covers 275 million ha worldwide (i.e., 20% of cultivated land) and accounts for 40% of world food production [1]. In the semi-arid area of the middle Niger valley (Niger) (Figure 1), irrigation techniques have been developed to respond to aridity and an increasing population [2]. Three factors have contributed to the rapid development of irrigation: (i) successive droughts from 1972 to 1973 and 1983 to 1984, which made people aware of serious risks to rainfed production; (ii) the high yields rapidly obtained in irrigated rice (*Oryza sativa*) and vegetable production; and (iii) commitment of the national government, farmers' organizations and several donors to irrigated agriculture. The first irrigated systems in Niger were constructed during the colonial period in the 1930s. From 1934 to 2011, 36 irrigated rice-growing systems were built along the middle Niger valley in Niger. Irrigation has improved the country's food security but has led to serious soil degradation by salinity. Irrigation facilities in many of these irrigated systems have aged and have not been renovated, and their drainage systems are not functioning. Excessive irrigation of less drained soil leads to waterlogging and soil salinization. The absence of a drainage network in clayey soils, which is characteristic of the study area [3], increases the concentration and precipitation of soluble salts at the soil surface but also in the subsoil and groundwater. Salt precipitates appear as white spots, mainly in the bright red oxidized clay horizons and in association with yellow iron spots. This salinization process threatens the sustainability of crop production in the study area by decreasing yields or reducing them to zero. The mechanisms by which soil salinity affects plant growth are generally known and have been summarized by many researchers [4–6]. In the middle Niger valley, rice is the main crop grown in irrigated fields. Rice can tolerate salinity without a reduction in yield up to a threshold of 3.0 dS/m of electrical conductivity of saturated paste (ECe) [7]. To preserve soil as a natural resource and maintain sustainable crop production in the study area, the spatial extent and level of salinity must be known. This can be done by performing field surveys to measure soil EC or electrical resistivity. Due to the large size of the Niger River basin (>85 km2), however, doing so would require large amounts of time, labor and money. Using the potential of remote sensing along with other data sources is currently a promising method to map salinity at a large scale with high accuracy.

**Figure 1.** Location of the study area and aerial view of the experimental fields in the Sébéri and Tchagriré irrigation systems south of Niamey (Niger).

The use of remote sensing to monitor and map soil salinity dates to the 1990s [8]. Many researchers have used remote-sensing and ground-truth data (i.e., direct or indirect measurements of soil) to assess and monitor salinity [9,10]. A variety of remote-sensing data have been used to identify and monitor salt-affected areas: aerial photographs, visible and infrared multispectral images, video images, infrared thermography, microwave images and data collected by airborne geophysics and electromagnetic induction [11]. Two remote-sensing methods exist to identify soil salinity: (i) observing surface conditions of bare soil (i.e., salt efflorescence) and (ii) monitoring the behavior of vegetation affected by salinity over time [8]. Before the most recent generation of satellites was developed, many researchers used a classic single-date approach [12] because most older satellites passed over a given location too infrequently to enable monitoring. However, this single-date approach has some limits: although highly saline and non-saline soils can be detected easily, intermediate salinity classes are difficult to differentiate. Moreover, if vegetation is present, bare soil cannot be seen, and vegetation behavior cannot be monitored with a single-date approach.

The most recent generation of satellites (Sentinel-2, Landsat 8, SPOT 6 & 7) offers high spatial resolution and frequent passes over a given location, which makes identifying salinity variations and monitoring vegetation behavior appear possible. A variety of remote-sensing indices, such as the salinity index (SI) [13], brightness index, normalized difference salinity index and normalized difference vegetation index (NDVI) [14], have been developed to estimate the salinity of bare soil and monitor vegetation behavior in saline environments. These indices have been combined with ground-truth data. For instance, soil salinity was measured by electromagnetic induction and then combined with multi-year Landsat 7 reflectance data to map salt-affected soils in the western San Joaquin Valley, California, USA [9]. Electrical conductivity (EC) ground-truth measurements and various Sentinel-2 spectral parameters were used to create more reliable soil salinity maps in the Ebinur Lake region, Xinjiang, China [15]. To monitor salt-affected areas in Turkey, researchers performed multi-temporal monitoring of salinity using field EC measurements and spectral indices derived from multi-year Landsat 5 and 8 satellite images [16]. In practice, most studies have focused more on detecting severely salt-affected areas than on detecting and monitoring slightly or moderately affected areas.

This study aimed to develop a method to detect potential areas of soil salinity using multi-spectral and high-resolution Sentinel-2 satellite images combined with field data using two complementary approaches: (i) observing salinity of bare soil using the SI and (ii) monitoring vegetation behavior from 2016–2019 (eight growing seasons) in the arid zone of the Niger River basin. The study involved a four-year time series of Sentinel-2 remote-sensing images and field measurements of biomass and topsoil characteristics of cultivated rice fields.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area corresponded to the irrigated systems of Sébéri and Tchagriré (6.5 km2), located in the Niger River basin 50 km southeast of Niamey, the capital of Niger (13◦16 35.32 N, 2◦21 31.81 E) (Figure 1). The climate of the area corresponds to the dry tropical zone of the Sudano–Sahelian type. Annual precipitation has high spatial, temporal and interannual variability and a general trend towards a southward shift of isohyets over the past 30 years. Mean annual precipitation is ca. 510 mm/year (standard deviation (SD) = 100 mm/year) (National Meteorological Service of Niger). Precipitation is irregularly distributed in space and time: it peaks in August (150 mm) and is lowest in October (22 mm) and May (20 mm). Mean monthly temperature is 36 ◦C during the hottest period (April), with a maximum temperature of 47 ◦C. Minimum mean monthly temperatures are 25 ◦C and are observed from December to January. Evaporation varies from 1700 to 2100 mm/year. The water deficit is thus large during the dry season and is accentuated by the Harmattan, a dry continental trade wind from the Sahara. Consequently, the study area has two seasons: a dry season from October to May and a wet season from June to September [17].

Paddy fields in Sébéri and Tchagriré are located on the lowest alluvial terrace of the Niger Valley (Figure 1). The soils are Vertisols (60–74% clay in topsoil and 52–85% in subsoil) but also acidic, with a pH in water of 4.0–6.2 (SD = 0.5), an EC of 0.01–7.20 dS.m−<sup>1</sup> (SD = 2.3 dS.m−1), and they have low hydraulic conductivity at water saturation, which ranges from 2.8 <sup>×</sup> 10−<sup>8</sup> in the topsoil horizon to 1.5 <sup>×</sup> <sup>10</sup>−<sup>8</sup> m s−<sup>1</sup> at a depth of 50 cm [18,19]. In Sébéri, EC generally increases from the areas next to the sand dunes to the river's floodwater protection dike [19]. The types of salts found are hexahydrite, gypsum, epsomite and, secondarily, wattevilleite and sodium carbonate hydrate [19].

The observed long-term downward trend in rice yields can be explained by several factors: poor management of agricultural equipment, lack of technical supervision, soil degradation due to salinity and failure to respect the cropping calendar. Implementing a fixed and common cropping calendar for all rice-growing systems would create two growing seasons per year during optimal climatic conditions. It provides for two harvests, one in the dry cropping season (mid-November to mid-May) and the other in the wet cropping season (mid-June to mid-December) (Figure 2).

**Figure 2.** Double-cropping calendar for rice in the middle Niger valley (adapted from [20]).

#### *2.2. Field Data Collection Strategy*

In the Niger River basin, rice is grown in 0.25 ha (25 m × 100 m) fields. We selected 44 and 20 fields in the Sébéri and Tchagriré irrigated systems, respectively, for data collection on two dates in the 2019 dry season, each of which corresponded to a date when the Sentinel-2 satellite passed over the study area. We chose the 64 fields based on the level of salinity in the study area and whether they were cultivated with rice or bare. Three sampling plots of 1 m<sup>2</sup> each were set up on one diagonal of each of the 0.25 ha fields selected to collect soil and phenological parameters.

The first data collection campaign was performed on 8 February 2019, near the start of the growing season, when fields were already flooded and covered by a layer of irrigation water a few centimeters thick. The phenology of rice plants was assessed in a quarter (0.25 m2) of each 1 m2 plot. Then, all aerial biomass of this quarter was cut, weighed and dried under laboratory conditions to determine dry matter. The topsoil (0–30 cm) was sampled in each plot using an auger, and one composite sample per field was obtained by carefully mixing the three elementary plot samples. Given the rapidity of auger sampling and the very low hydraulic conductivity of these clay soils, sampling under a water layer did not significantly change the EC and pH values of the soil, as shown in previous studies performed in these irrigated systems [19]. Consequently, salt losses during soil sampling are very low, ensuring good reliability of collected samples and good representativeness of soil salinity measurements. The composite samples were air-dried and ground to pass through a 2-mm sieve. Then, pH in water (pH1:5) and EC in water (EC1:5) were analyzed in the laboratory following ISO 10390 and ISO 11265, respectively.

The second data collection campaign was performed on 4 May 2019, during the harvest period, when fields were not flooded. Phenological assessment and soil sampling were identical to those

during the first campaign. Moreover, rice grain was collected from an undisturbed 0.25 m2 quarter to estimate grain yield. Soil and phenological parameters were analyzed statistically for each field.

#### *2.3. Remote-Sensing Data Collection*

Remote-sensing data and in situ observations were processed in multiple steps (Figure 3). First, optical images were obtained from Sentinel-2 satellites. The Sentinel-2 mission is a constellation of two satellites that are 180◦ out-of-phase on the same orbit, which increases the frequency of passes over a given location (i.e., every five days). Each satellite records 13 bands, with three spatial resolutions: blue, green, red and near-infrared bands at 10 m resolution; red-edge and mid-infrared bands at 20 m resolution; and aerosol and water-vapor bands at 60 m resolution. Sentinel-2 s utility thus lies in its high revisit frequency and high spatial resolution.

**Figure 3.** General data-processing flowchart.

All existing pre-processed Sentinel-2 images from January 2016 to December 2019 were downloaded. Pre-processing was performed by the MACCS (Multi-sensor Atmospheric Correction and Cloud Screening) processing chain, which consists of three successive steps: (i) cloud detection (using the satellite cirrus band), (ii) aerosol-thickness estimation and (iii) atmospheric correction. This processing chain, developed by CESBIO and CNES, provides ortho-rectified atmospheric-corrected surface reflectance images of 100 km × 100 km as a final product [21].

Overall, 157 preprocessed Sentinel-2 images were downloaded and then resampled at a spatial resolution of 10 m. Reflectance values that were missing due to cloud cover were estimated as the mean reflectance value of the two dates before and after the missing date. Finally, the study area of the irrigated systems of Sébéri and Tchagriré was extracted from the entire image.

A second step consisted of distinguishing the periods during which a plot was bare or covered with vegetation. To do this, the NDVI was first calculated per pixel (Table 1) [22], and then the mean and standard deviation of the NDVI of each plot's pixels was calculated. A pixel was considered to be part of a plot if its centroid fell within the plot's polygon. Since the plots measured 25 m × 100 m, each contained ca. 25 pixels. Based on the literature [23,24], plots with a mean NDVI ≤0.21 or >0.21 were considered as bare soil or vegetation, respectively.

**Table 1.** Characteristics and equations of the bare soil and vegetation indices selected. RED and NIR indicate red and near-infrared bands of the Sentinel-2 images.

Means and standard deviations of NDVI were used to create a 4-year time series of NDVI for each plot. For the periods when the plots were estimated to be bare, we calculated an SI [25] (Table 1). Otherwise, for the periods when vegetation was dominant, we created a time series of NDVI over the four years of study.

A few days before rice planting and during the first part of the vegetation development cycle, a water layer covers the soil. This layer is thin, no more than a few centimeters thick, and has high turbidity. The NDVI of this water layer ranges from 0.10–0.21 before planting and then decreases as the rice is transplanted and the vegetation canopy develops. In comparison, NDVI values of the river near the plots are negative and close to −0.3.

The last step consisted of calculating the time-integrated NDVI (TI-NDVI) [26,27] for each growing season for fields with vegetation. TI-NDVI for a given growing season was calculated as follows:

$$TI - NDVI = \sum\_{t=d1}^{d\epsilon} \left[ \frac{1}{n} \sum\_{t=1}^{n} \left( \left( \frac{NDVI\_{i,t} + NDVI\_{i,t+1}}{2} \right) - 0.21 \right) \times \left( d(t+1) - d(t) \right) \right] \tag{1}$$

where *d*1 and *de* are the day of the year of the start and end (harvest) of the growing season (the same for all fields), respectively; *d(t)* is the day of the year of date *t* for the set of days for which Sentinel-2 images are available during the growing season; *n* is the number of pixels within the field; and *NDVIi*,*<sup>t</sup>* is the NDVI of pixel *i* on date *t*. Units of TI-NDVI are expressed as NDVI.days.

Figure 4 illustrates the calculation of TI-NDVI for two successive growing seasons in 2019 derived from NDVI estimates based on existing Sentinel-2 images.

**Figure 4.** Example of normalized difference vegetation index (NDVI) dynamics for a given field for the dry and wet seasons of 2019 (derived from 42 cloud-free images among 60 existing Sentinel-2 images) and time-integrated NDVI (TI-NDVI) calculation for the two growing seasons.

#### *2.4. Multidimensional Analysis of the Data*

Principal component analysis (PCA) was performed using the FactoMiner package in R software [28] to analyze the multidimensional space of the remote-sensing and field data from 2019. Field data were pH and EC at the start and end of the growing seasons, aerial biomass and grain yields of the 64 fields. Qualitative variables were added: pH class, salinity class, land use and position in relation to the direction of river flow. Remote-sensing data were the mean SI of each field, TI-NDVI of each of the eight growing seasons and the number of seasons that each field was used to produce rice in dry and wet seasons. Remote-sensing data were considered as active variables and field data as illustrative variables.

PCA was followed by hierarchical clustering analysis (HCA) using the module HCPC (hierarchical clustering on principal components) [29] in R software [30] to create clusters of fields that behaved similarly. Finally, supervised classification for grids (SCG) was applied to 10-m resolution grids of the entire area that described the remote-sensing variables used in the PCA-HCA analysis (i.e., 8 grids with the TI-NDVI of each pixel for the 8 seasons and 1 grid with the SI of each pixel) to represent the spatial distribution of the clusters. Each pixel is assigned to the cluster with the shortest Mahalanobis distance, and the distance to the nearest cluster is analyzed to assess the quality of assignment to a cluster. SAGA-GIS 2.3.2 [31], implemented in QGIS 3.4.3, was used to perform the SCG with a maximum-likelihood algorithm.

#### **3. Results**

#### *3.1. NDVI Dynamics over the Eight Growing Seasons*

Figure 5 shows the dynamics of mean NDVI estimated from the 25 pixels within each of the three fields. The threshold of NDVI = 0.21 was selected to identify periods without live vegetation, and it enabled the identification of the short intercropping period after harvest when field irrigation was stopped, the soil plowed and irrigation started again before transplanting the next rice crop.

**Figure 5.** Example of dynamics of mean NDVI for three fields over eight growing seasons from January 2016 to December 2019: (**a**) Field CP3\_1, cultivated in season 2 but not in the other seasons; (**b**) Field CP11-1, always cultivated but with lower NDVI in dry seasons than wet seasons; and (**c**) Field TVP7, always cultivated and with high NDVI in both dry and wet seasons. For each field, mean NDVI and its 95% confidence interval were estimated from the field's ca. 25 pixels.

Non-cultivated fields had NDVI dynamics below the threshold or slightly above due to weed growth (Figure 5a). When cultivated during a growing season, some fields had lower NDVI due to constrained plant growth, which could have been due in part to soil salinity. Among the fields cultivated during both dry and wet seasons, two groups of cultivated fields were identified: (i) those with lower NDVI during the dry seasons than wet seasons (Figure 5b) and (ii) those with smaller differences between dry and wet seasons and generally high NDVI (Figure 5c).

#### *3.2. Spatial Variation in TI-NDVI over the Eight Growing Seasons*

For the 8 seasons and 44 irrigated fields monitored in Sébéri, TI-NDVI was always significantly lower during dry seasons than wet seasons (Figure 6a): mean TI-NDVI per season ranged from 15 to 19 NDVI.days during the four dry seasons and 29 to 34 NDVI.days during the four wet seasons. TI-NDVI also varied greatly among fields for a given season and showed systematic trends, with values frequently lower in some fields in the south and northwest of the irrigation system and significantly higher in fields in the center (Figure 6b). This spatial variation was lower in some seasons (e.g., seasons 2 and 6) but higher in others, e.g., season 7 in 2019, when field data were collected.

**Figure 6.** Temporal and spatial variations in the time-integrated normalized difference vegetation index (TI-NDVI) for the 44 fields of the Sébéri irrigation system: (**a**) mean and standard deviation of TI-NDVI for each season; (**b**) maps of TI-NDVI estimated for the eight growing seasons.

#### *3.3. Salinity Measured in the Field in 2019*

For the 64 fields monitored in 2019 during the dry season, EC1:5 ranged from 0.01 to 5.36 dS/m at the start of the season and 0.44 to 6.19 dS/m at the end of the season (SD = 1.16 and 1.33 dS/m, respectively) (Table 2). The pH ranged from 5.52 to 6.68 in January and 5.29–6.25 in June (SD = 0.44 and 0.52, respectively). Differences between means and SDs for the two dates were not significant for EC1:5 or pH. Total aerial biomass and grain yield at harvest varied greatly (Table 2), and 18 fields had no rice production.

**Table 2.** Statistics of field data for the dry growing season in 2019 (i.e., season 7) (n = 64, including 18 non-cultivated fields).

#### *3.4. Correlation between Remote-Sensing Data and Field Data during the 2019 Dry Season*

A strong and significant (*p* < 0.05) positive correlation was observed between TI-NDVI in season 7 and both total aerial biomass (r <sup>=</sup> 0.77, *<sup>p</sup>* <sup>&</sup>lt; <sup>2</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>14) and grain yield (r <sup>=</sup> 0.82, *<sup>p</sup>* <sup>&</sup>lt; 2.2 <sup>×</sup> <sup>10</sup><sup>−</sup>16) (Table 3), which confirms that TI-NDVI is a good indicator of rice vegetation growth and its final yield. Soil EC1:5 was negatively correlated with TI-NDVI at the start and end of the season (r = −0.38, *p* = 0.002 for the start, *p* = 0.002 for the end), with total biomass (r = −0.23, *p* = 0.062) and, significantly, with grain yield (r = −0.29, *p* = 0.019).

**Table 3.** Pearson correlations between the TI-NDVI indicator derived from Sentinel-2-images and soil (EC1:5, pH1:5) and crop indicators (total biomass, grain yield). Data were collected for 64 fields during the dry season of 2019 (i.e., season 7). SS and ES indicate the start and end of the growing season, respectively. Bold values indicate significant (*p* < 0.05) correlations.

#### *3.5. PCA and HCA Analysis of Remote-Sensing Data over the Eight Growing Seasons*

The first axis of the PCA was positively correlated with all of the TI-NDVI variables (Figure 7), which indicates that the spatial variability in TI-NDVI among fields was the main factor that influenced variations in the dimensional space. The second axis was correlated with SI but also distinguished the TI-NDVI estimates of the dry seasons (1, 3, 5 and 7) from those of the wet seasons (2, 4, 6 and 8).

While soil pH was little represented in the first PCA plane, soil EC1:5 measured at the start and end of season 7 was negatively correlated with the first axis and thus with TI-NDVI variables. More surprisingly, neither soil EC1:5 variable was correlated with SI. Following PCA, HCA identified five clusters of fields that behaved similarly during the eight growing seasons according to the remote-sensing data (Figure 8).

**Figure 7.** First plane of principal component analysis with variables derived from remote sensing as active variables (black) and soil and vegetation variables as illustrative variables (blue). Abbreviations: TI-NDVI\_SN: time-integrated normalized difference vegetation index for season N; Mean\_SI: mean salinity index; biom\_yield: biomass yield at harvest; grain yield: grain yield at harvest; Soil EC\_start\_s: soil electrical conductivity at the start of the season; soil EC\_end\_s: soil electrical conductivity at the end of the season; soil pH\_start\_s: soil pH at the start of the season; Soil pH\_end\_season: soil pH at the end of the season; Nbr\_dry\_s: number of dry growing seasons during the time series; and Nbr\_wt\_s: number of wet growing seasons during the time series.

**Figure 8.** Representation on the first plane of principal component analysis (PCA) of the five clusters of fields defined by ascending hierarchical classification applied to the first three axes of the PCA. Points are labeled with field codes starting with S or T respectively for Seberi or Tchagrire irrigation system, while squares indicate barycenters of clusters.

#### *3.6. Description of the Field Clusters*

The five clusters differed in their remote-sensing and field data (Table 4):

**Table 4.** Mean (and standard deviation) per cluster of fields defined by ascending hierarchical clustering for the variables derived from Sentinel-2 images (TI-NDVI, SI) and from field data collected in 2019 (EC1:5, pH, total biomass, grain yield). For each variable, different letters indicate significant differences in the mean among clusters according to a *p* < 0.05 Tukey test at a 95% confidence level.

### *3.7. Mapping the Clusters over the Entire Study Area*

Using the eight grids with the TI-NDVI of each pixel for the eight seasons and the grid with the SI of each pixel, SCG enabled each pixel to be attached to one of the five clusters defined from the 64 study fields considered as training areas. Figure 9 shows the distribution of the five clusters over the entire study area, with the associated distance to the nearest cluster, which indicates when the cluster represents a given pixel well (i.e., small distance) or poorly (i.e., large distance). Clusters 1 and 2, both of which had low TI-NDVI during the dry seasons since they were often not cultivated, were located mainly in the north and south of the Sébéri system (Figure 9). Cluster 5, with the highest TI-NDVI in both dry and wet seasons, predominated in the north of the Tchagriré system and the center of the Sébéri system. Cluster 3, with moderate TI-NDVI in both seasons, represented large areas in Sébéri in intermediate positions between clusters 1 or 2 and 5, but also in non-agricultural areas (e.g., paths, natural areas), which had large distances to the nearest cluster. Cluster 4 occupied small areas in Sébéri and the south of Tchagriré, often near fields of cluster 2.

**Figure 9.** Mapping of the five clusters over the entire study area of the irrigated systems of Sébéri and Tchagriré using supervised classification that assigned each pixel to the nearest cluster.

#### **4. Discussion**

This study continuously monitored spatio-temporal dynamics of SI and NDVI in paddy fields using Sentinel-2 satellite images over eight growing seasons from 2016 to 2019. Based on field data collected (biomass, EC and pH) during the 2019 dry growing season, a relation between spectral indices (NDVI and SI) and field data was established to understand the behavior of rice crops and relate the spatio-temporal variation and pattern of spectral indices to salinity.

#### *4.1. Variation in Spectral Indices among Crop Seasons*

In the 2016–2019 time series, the spectral indices used (SI and NDVI) varied in different ways. First, SI was estimated at dates when bare soils prevailed according to an NDVI threshold (NDVI <0.21). SI averaged over the 4 years of study varied in a narrow range and was weakly correlated with soil EC1:5; only highly saline areas, generally not cultivated, had higher SI values, due to the presence of salt crusts at the surface. These results may be explained by the short periods during which bare soil can be observed in irrigated paddy field systems; they generally last less than a month between successive crops, which corresponded to 3–5 dates when Sentinel-2 images were available. Thus, only a few dates when SI could be estimated were available. The soil may also be covered by crop residues or change drastically in water content due to the stopping or starting of irrigation after harvest or before preparatory work for planting rice. These factors may interfere with the estimation of SI and limit its ability to distinguish soil salinity.

In contrast, NDVI dynamics could be followed at a fine temporal resolution since cloud-free images were available at 133 dates during the eight growing seasons, with a mean of 21 dates per season since the end of 2017. We observed significant differences in TI-NDVI among the fields and eight growing seasons. TI-NDVI, which differed greatly between fields in a given season and over time, was used to differentiate areas with constraints to vegetation growth in the irrigated systems. In a given season, non-cultivated or cultivated areas with constraints had low TI-NDVI in dry and wet

seasons (clusters 1 and 2) (Table 4). Areas cultivated throughout the year with few constraints had high TI-NDVI in wet seasons and moderate TI-NDVI in dry seasons (cluster 3), while zones cultivated in dry and wet seasons without any constraints had high TI-NDVI in both seasons (cluster 5) (Figure 8). Nonetheless, there were areas near the main drainage channel with salinity constraints and moderate TI-NDVI that were cultivated only in wet growing seasons. They may have occurred because they are always wet and wild grasses grow there in the dry seasons, which influences the NDVI. TI-NDVI also differed significantly between dry and wet seasons (Figure 5a,b). Wet growing seasons had higher mean TI-NDVI than dry ones. TI-NDVI varied from 29 to 34 NDVI.days in wet seasons and 15 to 19 NDVI.days in dry seasons. This result is due to the climatic conditions in the study area in dry seasons (hot, dry wind, high temperature) that limit crop growth. In addition, high temperature favors the rise of salt from lower horizons to the rooting zone of crops, which damages their roots. In wet seasons, conditions are more favorable for crop growth. Over the time series, wet and dry seasons varied among years, again due to climatic conditions (temperature in dry seasons and precipitation in wet seasons).

#### *4.2. Temporal and Spatial Patterns of NDVI*

Soil salinity influences vegetation density and crop growth, as explained by NDVI. TI-NDVI was used in this study to indicate constraints to crop growth and density, as in previous studies [32]. Five clusters were obtained to differentiate areas with constraints to rice growth in the two irrigated systems. Non-cultivated areas had the lowest TI-NDVI and were considered likely to be areas with high salinity constraints (cluster 1) (Figure 9). Areas cultivated only in wet seasons because of constraints had low TI-NDVI and may be considered to have moderate salinity constraints that limit crop growth (clusters 2 and 4) (Figure 9). Areas cultivated in both seasons with few constraints had high TI-NDVI and may be considered to have few salinity constraints (cluster 3) (Figure 9). Finally, areas cultivated without constraints had the highest TI-NDVI (cluster 5) (Figure 9) and were considered zones without salinity constraints. Our results show that the spatial pattern of TI-NDVI corresponds to the spatial distribution of problematic patches in the study area. The main constraint may be salinity, while other constraints could be explained by soil type, microtopography and farming practices. Mean ECe measured during ground truthing was calculated for each cluster and assigned to the respective TI-NDVI of the clusters to validate the four salinity classes obtained from the classification of TI-NDVI. Based on EC measurements, visual observations and knowledge of the terrain [33], mean TI-NDVI of the clusters were classified into four classes of soil salinity [34]: non-saline, slightly, moderately and very saline (i.e., ECe < 2, 2–4, 4–8 and 8–16 dS/m, respectively). Measured soil EC alone could not explain the level of salinity of the clusters (Table 4); however, integrated interpretation using the EC, mean TI-NDVI, total biomass and grain yield with visual observations of the study area could be used. Cluster 1, consisting of abandoned fields with low mean TI-NDVI in dry and wet seasons, no rice biomass or grain yield, maximum EC of 5.36 dS/m<sup>2</sup> and white efflorescence on the soil surface, can be considered saline soil. Clusters 2 and 4, consisting of fields cultivated only in wet seasons and rarely in dry seasons, despite having a lower level of surface salinity, have low TI-NDVI in dry seasons, low biomass and grain yield in the few fields cultivated in dry seasons and no yield in non-cultivated fields. Salt efflorescence is present in these fields. The maximum EC was 2.59 dS/m2, which may not reflect the true EC of the clusters. These two clusters can be considered moderately saline. Clusters 3 and 5 consisted of cultivated fields in both dry and wet seasons. Based on measured EC, they were classified as non-saline, but since cluster 3 has lower TI-NDVI, biomass and grain yields than cluster 5, cluster 3 has few salinity constraints. Thus, cluster 3 can be considered slightly saline and cluster 5 non-saline.

#### *4.3. Field EC Variation and Salinity*

In the Sébéri irrigated system, EC generally increased from the areas next to the sand dunes to the river's floodwater protection dike, perhaps because the dike has modified the functioning of the soils next to it [18]. In the Tchagriré irrigated system, EC generally increased from the areas next to the river to the main drainage ditch, likely due to topography. Areas at lower elevations had higher EC than those at higher elevations. Overall, the pattern of EC followed those of TI-NDVI, biomass and grain yield. EC measured in the two systems does not reflect the real situation in the field. Areas that showed signs of salinity (i.e., abandoned fields, low NDVI, low yields) had low measured EC, which indicates that a salt stock may lie in the lower horizons. To map salinity in this situation, measured EC should be compared to other maps (e.g., TI-NDVI, yield, soil, elevation) of the same site with similar sampling patterns or resolution. Doing so may provide useful insights into other parameters that could explain salinity.

#### **5. Conclusions**

This study developed a step-by-step method to estimate constraints on the growth of rice, the most important of which is salinity. Dense time series of Sentinel-2 images over eight growing seasons enabled us to describe the behavior of rice biomass and to differentiate fields where rice is subjected to stress during growing seasons. In irrigated systems, periods of bare soil were brief, and the SI derived from Sentinel-2 images could not differentiate soil salinity of fields. Monitoring vegetation behavior over four years by deriving the NDVI from Sentinel-2 images and calculating the TI-NDVI was able to differentiate fields based on constraints that limit rice growth. Several constraints can occur in areas subjected to stresses but can be related locally to soil salinity and verified by field sampling, which can be guided by TI-NDVI classification. This approach is particularly adapted to irrigated rice systems in which monoculture prevails and differences in TI-NDVI are not caused by different crops.

**Author Contributions:** I.M., C.W. and D.M. conceived the topic and developed the method. P.P. performed the technical work of downloading and processing images and calculating spectral indices. I.M., C.W. and D.M. wrote the manuscript, which was improved by H.N., Y.G. and I.A.B. participated actively in the fieldwork. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was made possible thanks to financial support of the EU-funded Agrocampus ERASMUS+ cooperation project in Rennes, France. This program financed a one-year's internship for M.I. and his supervision.

**Acknowledgments:** This study was performed with the administrative and technical support of the soil science laboratory from UMR SAS, INRAE and Institut Agro.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Framework for Accounting Reference Levels for REDD+ in Tropical Forests: Case Study from Xishuangbanna, China**

**Guifang Liu 1,2, Yafei Feng 1, Menglin Xia 1, Heli Lu 1,2,3,\*, Ruimin Guan 1, Kazuhiro Harada <sup>4</sup> and Chuanrong Zhang <sup>5</sup>**

**Abstract:** The United Nations' expanded program for Reducing Emissions from Deforestation and Forest Degradation (REDD+) aims to mobilize capital from developed countries in order to reduce emissions from these sources while enhancing the removal of greenhouse gases (GHGs) by forests. To achieve this goal, an agreement between the Parties on reference levels (RLs) is critical. RLs have profound implications for the effectiveness of the program, its cost efficiency, and the distribution of REDD+ financing among countries. In this paper, we introduce a methodological framework for setting RLs for REDD+ applications in tropical forests in Xishuangbanna, China, by coupling the Good Practice Guidance on Land Use, Land Use Change, and Forestry of the Intergovernmental Panel on Climate Change and land use scenario modeling. We used two methods to verify the accuracy for the reliability of land classification. Firstly the accuracy reached 84.43%, 85.35%, and 82.68% in 1990, 2000, and 2010, respectively, based on high spatial resolution image by building a hybrid matrix. Then especially, the 2010 Globeland30 data was used as the standard to verify the forest land accuracy and the extraction accuracy reached 86.92% and 83.66% for area and location, respectively. Based on the historical land use maps, we identified that rubber plantations are the main contributor to forest loss in the region. Furthermore, in the business-as-usual scenario for the RLs, Xishuangbanna will lose 158,535 ha (158,535 <sup>×</sup> <sup>10</sup><sup>4</sup> m2) of forest area in next 20 years, resulting in approximately 0.23 million t (0.23 <sup>×</sup> <sup>10</sup><sup>9</sup> kg) CO2e emissions per year. Our framework can potentially increase the effectiveness of the REDD+ program in Xishuangbanna by accounting for a wider range of forest-controlled GHGs.

**Keywords:** reference levels; REDD+; greenhouse gas emissions; Xishuangbanna; monitoring and reporting

#### **1. Introduction**

Forests account for almost half of the global terrestrial carbon pool, and the vegetation within them alone (excluding soils) holds approximately 75% of all living carbon. The total carbon content in forest ecosystems is estimated to be 638 Gt [1–5]. Tropical forests play a particularly important role in the global carbon budget because they contain as much carbon in their vegetation and soils as all the temperate-zone and boreal forests combined [6–12]. Per unit area, tropical forests store, on average, approximately 50%

**Citation:** Liu, G.; Feng, Y.; Xia, M.; Lu, H.; Guan, R.; Harada, K.; Zhang, C. Framework for Accounting Reference Levels for REDD+ in Tropical Forests: Case Study from Xishuangbanna, China. *Remote Sens.* **2021**, *13*, 416. https://doi.org/ 10.3390/rs13030416

Received: 30 December 2020 Accepted: 20 January 2021 Published: 26 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

more carbon than their nontropical counterparts. Scientists agree that to achieve the goals of the United Nation's Framework Convention on Climate Change (UNFCCC), namely avoiding irreversible damage to the climate system, global warming must not exceed 2 ◦C [13–16]. However, concentrations of CO2 in the atmosphere are already so high that global emissions will likely peak before they start to decline. Thus, in order to remain under the above-mentioned threshold, emissions from all major sources (i.e., from developed countries, major developing country emitters, and deforestation) must begin to decline within the next decade [17–19].

The Conference of the Parties (COP) agreed that Reducing Emissions from Deforestation and Forest Degradation (REDD+) with the enhancement of the removal of greenhouse gas (GHG) emissions by forests in developing countries could support the goals of the framework through the positive incentives provided by the UNFCCC. The general consensus at Doha in 2012 following last year's COP17 was that the results of financing, safeguards, measurements, and reporting and verification for REDD+ were mixed. In addition, significant progress has been widely recognized as having been made only within the technical arena relating to reference levels (RLs). The 19th Conference of the Parties to the UNFCCC (COP19) and the 9th Conference of the Parties to the Kyoto Protocol (CMP9) were jointly held on the topic of REDD+ funding in Warsaw, Poland, and in-depth discussions on the action points were conducted.

GHG-based compensation for REDD+ requires an agreement on emission RLs. Key elements for setting these RLs include the ability to measure changes throughout all forested areas, the use of consistent methodologies at repeated intervals to obtain accurate results, and the verification of results with ground-based or very high-resolution observations [20–24]. RLs have profound implications for the effectiveness of climate-related policies, cost efficiency, and distribution of REDD+ financing, and they involve a number of tradeoffs [25–30]. In this paper, referring to the business-as-usual scenario, we introduce a methodological framework for setting RLs for REDD+ applications in tropical forests in Xishuangbanna, China, by coupling the Good Practice Guidance (GPG) on Land Use, Land Use Change, and Forestry published by the Intergovernmental Panel on Climate Change (IPCC) and land use scenario modeling. This study contributes to the literature by highlighting key challenges for setting RLs as part of the REDD+ program.

#### **2. Data and Methodology**

#### *2.1. Research Area*

Not only is the forest in the Xishuangbanna region (Figure 1) the world's largest preserved area located in the northernmost part of the Earth, but it is home to the majority of tropical forest ecosystems in China as well. The geology, climate, and soil of Xishuangbanna are suitable for the growth and reproduction of various organisms. Moreover, 4500 species of higher plants have been recorded in Xishuangbanna, accounting for about one-seventh of the total number of higher plants in China. The native vegetation types include those found in tropical rain forests, montane rain forests, tropical monsoon forests, subtropical evergreen broad-leaved forests, deciduous broad-leaved forests, warm coniferous forests, and bamboo forests as well as shrubs and grasses [31–34]. In recent years, due to the increase in the population, intensification of anthropogenic activities, the enabling climate, and suitable terrain conditions in the area, the cultivation of rubber and other tropical and economically important crops has risen rapidly. Thus, the changes to the forest have been very dramatic.

**Figure 1.** Landscape of the Xishuangbanna, China.

#### *2.2. IPCC's Good Practice Guidance*

The IPCC's existing GPG for Land Use, Land Use Change, and Forestry provides the recommended approach to account for fluctuations in carbon stocks resulting from changes in the use and management of forests. This framework has been accepted by all Parties in the Bali Action Plan of COP13 [35–37]. The IPCC's GPG framework refers to two basic inputs for forest carbon accounting, namely activity data and emission factors. Activity data in the REDD+ context refer to the areal extent of emissions. For example, in the context of deforestation, activity data refer to the area of deforestation, presented in hectares (104 m2) over a known time period. Emission factors refer to the emission or removal of GHGs per unit activity. The emission or removal of GHGs resulting from land use conversion ultimately alters ecosystem carbon stocks.

#### 2.2.1. Emissions Factors

To estimate emissions factors, the required number of sample plots was determined to the necessary accuracy using the size of the forest area and other available resources. Provisional surveys and/or existing data can be utilized to establish sample sizes, and tools also exist to calculate sample sizes based on fixed precision levels or given fixed inventory costs [38–41]. In the event carbon stocks and flows are to be monitored over the long term, permanent sites should be considered in order to reduce between-site variability and to capture actual trends as opposed to short-term fluctuations [42].

In the study, the aboveground biomass density map was sourced from Global Forest Watch (http://www.globalforestwatch.org/). This map is a global aboveground biomass density map produced in 2000 according to the method devised by Baccini [43]. Based on the improved methodology, the resolution can be increased to as much as 30 m. The aboveground biomass density map of Xishuangbanna region was extracted using the mask extraction method.

More recently, Maurizio Santoro [44] have proposed an integration methodology for estimation of aboveground biomass density for around the year 2010 by combining SAR, LiDAR, and optical observations together with other datasets such as auxiliary datasets

from forest inventories, additional remote sensing observations, climate variables, and ecosystems classifications. We compared with the latest biomass map developed by Santoro against Global Forest Watch product.

#### 2.2.2. Activity Data

Estimation of activities associated with national-level deforestation monitoring is practically possible only via remote sensing [45–49]. Since the early 1990s, changes in forest area have been monitored from space with confidence. Some countries have had well-established operational systems for over a decade.

Taking into account the availability of the data and their matching, this study used the Enhanced Thematic Mapper Plus/Thematic Mapper (ETM+/TM) remote sensing images related to Path number 130 and Row number 044, Path number 131 and Row number 045, Path number 130 and Row number 045, Path number 129 and Row number 045 in the study area in 1990, 2000, and 2010 to interpret land use changes(Table 1). The TM/ETM+ remote sensing images and normalized difference vegetation index (NDVI) data were sourced from the US Geological Survey (USGS, http://earthexplorer.usgs.gov/). We mosaiced the TM/ETM+ remote sensing images of the same year, used ENVI to perform geometric correction and radiometric correction, converted all the map data projections to WGS84/UTM Zone47N (EPSG: 32647), used the Xishuangbanna administrative vector map for mask extraction, and performed cropping to obtain the images of Xishuangbanna.

**Table 1.** Landsat imagery used in this study.

The vegetation in the study area changes obviously with the seasons, and the NDVI values at different times, thus, have a greater influence on the research results. Thus, this study synthesized the maximum value of the NDVI data in multiple phases of the same year and also eliminated the influence of cloud cover on the research results.

According to the characteristics of land use cover in the study area (Table 2), we divided the land use cover into eight types: forestland, shrubland, grassland, cultivated land, rubber forest, tea gardens, construction land, and water. The training samples were determined using QuickBird images in Google Earth. The terrain of the study area is relatively complex, and many "homogeneous spectrum" phenomena occur in the interpretation of remote sensing images. To avoid this phenomenon, it is necessary to select as many training samples as possible. Different band combinations of Landsat7 ETM images have different characteristics. The selection of the training samples was carried out according to these characteristics. The ETM541 band combination is helpful for distinguishing different vegetation types when supplemented by NDVI data. The training samples of natural forests, shrubs, rubber plantations, and tea gardens were selected. The ETM453 band combination was used to select the cultivated land and water, while the construction land was extracted through the ETM743 band combination, and the remainder was categorized as other land. Supervised classification was performed using the selected training samples to obtain the preliminary classification results, and the accuracy test was conducted. If the results did not agree, the training samples were reselected, and the supervised classification and accuracy tests were reperformed until they were ideal. Finally, the classification results were recoded, clustered, and eliminated, and the broken patches were merged into the adjacent largest classification to unify the smallest unit.

**Table 2.** Land use classification followed for this study.

The Land Use Dynamic Index considers the transfer of land use types during the study period, and reflects the intensity of regional land use changes during this time. It is essential to find hot spots of land use changes at different spatial scales. It is one of the important parameters to analyze the dynamic changes in land use space [50]. Equation (1) was used to calculate the index.

$$K\_i = \frac{\mathcal{U}\_{t\_1} - \mathcal{U}\_{t\_2}}{\mathcal{U}\_{t\_1}} \times \frac{1}{t\_2 - t\_1} \times 100\% \tag{1}$$

*Ki* is the land use dynamic degree for land use type *i* in a certain period of time, *Ut*<sup>1</sup> and *Ut*<sup>2</sup> are the number of certain land use types at the start of the period *t*<sup>1</sup> and its end *t*2, respectively, and *t*<sup>2</sup> − *t*<sup>1</sup> is the research duration.

$$S = \left[ \sum\_{i=1}^{n} \left( \frac{\Delta S\_{i-j}}{S\_i} \right) \right] \times 100 \times \frac{1}{t} \times 100\% \tag{2}$$

*S* is the comprehensive land use dynamic degree in the study area corresponding to *t* time period. Δ*Si*−*<sup>j</sup>* is the area of land use type converted *i* converted to other land use types in the study period; *Si* is the area of type *i* land use type at the beginning of the study; *t* is the time period of land use change.

#### *2.3. Land Use Scenario Modeling for Reference Levels*

Land use simulation is based on years of known land use changes. It predicts future land use changes. Most land use models used to simulate the process of land use change typically need to solve two problems: the quantity problem and the distribution problem. The quantity problem refers to how much of the land area has changed, while the distribution problem involves pinpointing where those land changes occurred. This study applied the Land Change Modeler (LCM) [51–53], which uses the Markov chain model to predict the number of future land use changes, and then calculates the distribution location of these changes according to the Multilayer Perceptron (MLP) model.

Markov chain is a kind of "no after-effect" random stored procedure, as it assumes that the state of the current variable is only related to its previous state, not to its states at other moments. Therefore, it has good operability and is used in the simulation of various land use changes. In Equation (3) of Markov chain, for any positive integer *n* and possible states *i*0, *i*1, ..., *in* of the random variables,

$$P(X\_n = i\_n | X\_{n-1} = i\_{n-1}) = P(X\_n = i\_n | X\_0 = i\_0, X\_1 = i\_1, \dots, \dots, X\_{n-1} = i\_{n-1}) \tag{3}$$

As the land use change conforms to the basic characteristics of the Markov process, it can be regarded as a Markov process. Therefore, the Markov chain analysis can describe the land use change process and predict the future land use change trend. It is an important transformation tool in land use change modeling. However, the following prerequisites must be fulfilled [54–56]: (1) In a certain area, different types of land use should be transformable into each other, (2) the conversion between different types of land use can include many events, which are difficult to describe with a specific formula, and (3) within

the time limit of the study, the conversion status of the land use structure is relatively stable, which meets the requirements of the Markov chain. Moreover, the area ratio of the mutual conversion between the types of land uses equals the state transition probability.

MLP is a very widely used neural network in remote-sensing image processing, especially remote sensing image classification. MLP was used in the model primarily to calculate the land use change potential, that is, the future conversion probability between each land use type. The process involved analyzing future land use by establishing a land use driving force model and the quantitative relationship between each land use type to assess the probability of change. Based on the calculated potential distribution of soil use changes, the location of possible future land use changes can be determined. The back propagation algorithm used in MLP consists of two parts, namely the forward propagation of information and the backward propagation of errors. In the forward propagation process, the input information is calculated from the input layer through the hidden layer to the output layer, and each layer for the state of a neuron only affects the state of the next layer of neurons [57,58]. If the expected output is not obtained in the output layer, the error change value of the output layer is calculated, and then turned to reverse propagation, and the error signal is returned back along the original connection path through the network to modify the weights of neurons in each layer until the desired target is reached. During the forward propagation process, the state of the activated neuron is updated layer-by-layer from the input layer to the output layer, as shown in Equation (4):

$$\alpha\_{\vec{j}} = \sum\_{i} a\_{i} w\_{\vec{j}i}.\tag{4}$$

*xj* represents the total input received by neuron *j*, *wji* represents the weight between neurons *j* and *i*, *α<sup>i</sup>* denotes neuron *i* once *xj* is calculated. The most commonly used mapping function is the S (sigmoid) function, as shown in Equation (5).

$$a\_{\vec{\jmath}} = \mathbf{f}(x\_{\vec{\jmath}}) = \frac{1}{1 + \frac{1}{\exp\left(\frac{x\_{\vec{\jmath}}}{\vec{\jmath}}\right)}}.\tag{5}$$

It is crucial to check the accuracy and effect of the model to determine whether the model needs to be adjusted. The Receiver Operating Characteristic (ROC) curve test evaluates the model by comparing the predicted land change probability distribution map with the actual changed 0–1 map (the changed land value is 1, and the unchanged land value is 0) [59–61]. This step converts the simulated and reference images into a 2 × 2 table, with each table corresponding to a different threshold. The number of pixels within the thresholds of A, B, C, and D create the statistical figure for each ROC curve threshold. The following data are produced; *x* and *y* form the point (*x*, *y*), where *x* is the ratio of classifications labeled as true−, namely D/(B + D), and *y* is the proportion of true+ classifications, that is, A/(A + C). In order to be expressed as a positive value on the *x* axis, the opposite part of true− is generally represented by B/(B + D). Thus, the ROC curve test provides the Area Under the ROC Curve (AUC), which is obtained using the following formula:

$$\text{AUC} = \sum\_{i=1}^{n} (x\_i - x\_{i+1}) \times \left\{ y\_i + \frac{y\_{i+1} - y\_i}{2} \right\} \tag{6}$$

where *xi* refers to *x* for each threshold *i*, that is, B/(B + D), and *y* is calculated using D/(B + D).

#### **3. Results and Discussion**

#### *3.1. Analysis of Historical Land Use*

The land use maps in 1990, 2000, and 2010 and the accuracies are shown in Figure 2 and Appendix A. Firstly we randomly generated 2866, 2549, and 2481 sample points in 1990, 2000, and 2010 through hierarchical random sampling method. There were 1520 sample

points, 1227 sample points, and 1008 sample points in forest area in 1990, 2000, and 2010, respectively. Then we evaluated the accuracy of classification for sample points based on Google earth.

**Figure 2.** Land use maps and the accuracies in 1990, 2000, and 2010.

GlobeLand30, which was developed by the National Geomatics Center of China, is an open-access 30m resolution global land cover data product with an overall classification accuracy of over 80% [62,63]. We compared the area and spatial location of the forest land in 2010 extracted by Globeland30 with those in this study (Figure 3). Firstly, about 600 sample points are randomly generated within Xishuangbanna administrative region. Then these sample points are overlapped with Globeland30 and land use map respectively. Finally we evaluate the accuracy of land use map based on the consistency of forest land and nonforest land in Globeland30.

**Figure 3.** (**a**) Forest area in 2010 from Globeland30 and this study; (**b**)spatial distribution of the validation samples in Globeland30; (**c**) spatial distribution of the validation samples in this study.

In terms of the forest land area, the forest land area extracted from Globeland30 was 1.21× <sup>10</sup><sup>6</sup> ha, and that from this study is 1.05×106 ha, with the accuracy 86.92%. In terms of spatial location of the forest land, among 600 randomly generated sample points, 376 were the forest land and 230 were nonforest land in Globeland30; in comparison, 331 sample points were the forest land and 275 sample points were nonforest land in this study. The overall accuracy is 83.66%, and the kappa coefficient is 0.657.

The areas, changes, and dynamics of the three types of land use in 1990, 2000, and 2010 are shown in Figure 4.

**Figure 4.** *Cont*.

**Figure 4.** (**a**)Land use areas in 1990, 2000, and 2010, (**b**) changes in land use type, and (**c**) land use change degree in Xishuangbanna between 1990 to 2000, 2000 to 2010, and 1990 to 2010.

As shown in Figure 4, the areas under cultivated land, forested land, and water bodies in Xishuangbanna showed a downward trend from 1990 to 2010. Among them, the decrease in forest area is the most obvious, with a total reduction of 360,819 ha (360,819 × 104 m2) over the past 20 years, a dynamic land use change degree of −1.42%, a decrease of 265,491 ha (265,491 × 104 m2) from 1990 to 2000, and a reduction of 95,328 ha (95,328 × 104 m2) from 2000 to 2010. The area of cultivated land showed an increasing trend in the previous 10 years, marked by a rise of 18,153 ha (18,153 × 104 m2) and a dynamic land use change degree of 1.24%. The area of cultivated land decreased by a total of 21,456 ha (21,456 × <sup>10</sup><sup>4</sup> m2) in the latter 10 years, with a dynamic land use change degree of −1.31%. The area under water bodies declined continuously for the two decades, with a total reduction of 3996 ha (3996 × <sup>10</sup><sup>4</sup> m2) and a dynamic degree of −2.17%. Grasslands, rubber plantations, shrubland, tea gardens, and construction land in Xishuangbanna region showed increasing trends from 1990 to 2010. Among them, the area of rubber plantations showed the most obvious growth, with a total increase of 249,948 ha (249,948 × 104 m2) in 20 years, and a dynamic land use change degree of 9.87%. Moreover, the area under tea gardens increased by 43,686 ha (43,686 × <sup>10</sup><sup>4</sup> m2) in the past 20 years, the dynamic land use change degree being 6.26%. Although the areas under grassland, shrubland, and construction land increased, the changes were relatively insignificant.

In summary, the economic development of the Xishuangbanna region and the improvement in people's quality of life led to a rise in the cultivation of cash crops such as rubber and tea in the region in the past 20 years, resulting in a large number of forests being felled.

During the period 1990–2000, carbon emissions for Global Forest Watch and Santoro datasets were 7.85 million t CO2e and 5.63 million t CO2e, respectively, with a difference of 28.30%. During the period 2000–2010, carbon emissions for Global Forest watch and Santoro datasets were 2.82 million t CO2e and 2.00 million t CO2e, respectively, with a difference of 28.81%. Carbon emissions for the period 1990–2000 were about 2.8 times as much as those for the period 2000–2010 (Figure 5).

**Figure 5.** Carbon emissions for the period 1990–2000 and the period 2000–2010.

#### *3.2. Influencing Factors of Land Use Change*

There are many drivers that lead to deforestation and forest degradation within REDD+. Direct drivers are human activities or immediate actions that directly impact forest cover and loss of carbon such as agriculture expansion (both commercial and subsistence), infrastructure extension, and wood extraction. Indirect drivers are complex interactions of social, economic, political, cultural, and technological processes to cause deforestation or forest degradation. They act at multiple scales: international (markets, commodity prices), national (population growth, domestic markets, national policies, governance), and local circumstances (subsistence, poverty) [64–67]. Since RLs refer to the business-as-usual scenario, which means without any change in REDD+ drivers (situation, government, socio-economic forces, etc. that occur over time), this study only considered seven factors influencing land use change, namely distance to a road, distance to a river, elevation, slope, aspect, distance to an administrative center, and nature reserves (Table 3 and Figure 6).

**Table 3.** Factors influencing land use change and data acquisition methods.

**Figure 6.** Influencing factors of land use change in Xishuangbanna region.

Cramer's V coefficients (Table 4) were calculated to measure the correlation between the above-mentioned factors impacting land use change and land distribution. The larger the value, the stronger the correlation.

**Table 4.** Cramer's V coefficients indicating correlations between the influencing factors of land use change and land distribution.

#### 3.2.1. Distance to a Road

Besides playing a very important role in the economic and social development of a region, traffic conditions impact the land use status of a region. The overall correlation between the land type and distance from a road is 0.1334. Firstly, compared with the overall value, Cramer's V coefficient for shrubland and tea gardens is 0.2735 and 0.2148, respectively, which is much higher than the overall value. Thus, the distance from a road is a relative important factor affecting shrubland and tea gardens. Secondly, Cramer's V coefficient of the impact of the distance from a road on rubber plantations and construction land is 0.1550 and 0.1411, respectively, quite similar to the overall value. Thus, the affected land types dominated by road traffic in the Xishuangbanna region are shrubland, tea gardens, rubber plantations, and construction. It is evident that these land types are affected by anthropogenic activity. The reason of highest correlation between the shrubland and road is that it is very common in Xishuanbbanna to have roads built across shrubland rather other areas.

#### 3.2.2. Distance to a River

The precipitation in Xishuangbanna region is abundant and evenly distributed. The dependence of most land use types on rivers is not obvious, except for shrubland and tea gardens. Among them, the influencing factor, namely the overall correlation value of the distance from a river to the land type is 0.0905, and the Cramer's V coefficients for tea gardens (0.1277) are higher than this overall value. This result indicates that the distance from a river is the main factor affecting tea gardens.

#### 3.2.3. Terrain-Related Factors

Topographic factors play a very important limiting role in various production activities. The study area is mainly mountainous, and, thus, the topographic factors of elevation, slope, and aspect cannot be ignored. Firstly, the overall value of the correlation is 0.2539, and woodland and rubber plantations alone show higher correlation coefficients than this overall value (the corresponding Cramer's V coefficients are 0.3482 and 0.5297). During the period 1990–2010, the rubber plantation in Xishuangbanna continuously expanded from low-altitude flat valleys to mountainous areas in high altitudes due to high rubber price from the international market, population pressure, and economic development. This is the reason for the highest correlation between the elevation and rubber plantation. Cramer's V coefficients of elevation for shrubland, grassland, cultivated land, tea gardens, construction land, and other land are 0.1192, 0.1289, 0.1355, 0.1230, 0.0451, and 0.2130, respectively, indicating that their correlation coefficients are lower than the overall value.

Secondly, the slope affects the water distribution, wind speed, and soil texture required for crop growth. The overall value of the correlation for the slope is 0.1608, while Cramer's V coefficients for shrubland, rubber plantations, and tea gardens are 0.2131, 0.1984, and 0.1665, respectively, higher than the overall value. Thus, this factor can be regarded as the main factor impacting these land uses. However, in overall terms, Cramer's V coefficient is less than the corresponding values for grassland, cultivated land, construction land, and other land (0.0406, 0.0210, 0.1076, and 0.0256, respectively).

Finally, the aspect primarily affects the length of time and temperature for the growth and final yield of crops. The overall value in this case is 0.0431, while Cramer's V coefficients for shrubland, grassland, and rubber plantations are all greater than the overall value (0.0729, 0.0560, and 0.0570, respectively).

#### 3.2.4. Distance to an Administrative Center

Governmental administrative organizations are typically located in townships. Given the increasingly strict forest protection policies being applied to Xishuangbanna region, areas closer to governmental administrative organizations can be conveniently supervised and regulated, resulting in a certain deterrent effect on forest destruction and illegal mining of local resources. The overall value of the distance from a township is 0.1252. The corresponding Cramer's V coefficients for rubber plantations, and tea gardens (0.1542, and 0.1516, respectively) are higher than the overall value. However, the coefficients for grassland, cultivated land, construction land, and other land (namely, 0.0933, 0.0455, 0.1025, and 0.0361, respectively) are less than the overall value. Therefore, rubber plantations and tea gardens are clearly (and expectedly) impacted by distance to a township, whereas this is not so for the remaining land use types.

#### 3.2.5. Limiting Factor (Nature Reserve)

Xishuangbanna Nature Reserve is a national nature reserve consisting of five small subreserves, namely the Mengyang, Menglun, Mengla, Shangyong, and Manzhang Reserves. These sub-reserves are not geographically connected to each other and cover a total area of 242,500 ha (242,500 × <sup>10</sup><sup>4</sup> m2). Notably, 12.68% of the total area of the Prefecture is allocated to nature conservation, namely the protection of the tropical forest ecosystem and its rare wildlife. Relatively little land change has been observed in the protected area, and man-made damage has also been effectively contained. In this study, the conversion rate of certain land use types, such as forestland, in the protected area was set to 0; in other words, anthropogenic activities in these areas are completely restricted.

#### *3.3. Future Land Use Simulation Results and Inspection*

The expansion of rubber and other cash crops has caused massive forest loss and fragmentation in Xishuangbanna. The region experienced the most severe forest losses and degradation particularly for the period 1990 to 2010. Therefore, we chosen the period 1990 to 2010 for REDD+ in Xishuangbanna as the baseline, which is crucial to measure the emission reduction performance and consequently to negotiate meaningful deforestation emission reduction targets. As a result, the land use change data for 1990 and 2000 were used as inputs to the model of the Markov chain and MLP, and the 2010 land use change data were used as the verification values to simulate future land use. The validation of AUC value from the ROC curve method is 0.8, indicating that the results provided by the model are ideal. The land use prediction results for the Xishuangbanna region in the next 20 years of 2016–2035 are shown in Figure 7.

Area under forestland shows a downward trend and is the largest change over the 20 years, with the areal reduction amounting to 158,535 ha (158,535 × 104 m2). Conversely, the areas under rubber plantations, tea gardens, and cultivated land increase, with rubber plantations showing the highest increase (by 108,450 ha (108,450 × <sup>10</sup><sup>4</sup> m2)). The areas under tea gardens and cultivated land also increase, but only slightly (by 39,204 ha (39,204 × <sup>10</sup><sup>4</sup> <sup>m</sup>2) and 31,707 ha (31,707 × 104 <sup>m</sup>2), respectively). The areas under shrubland, grassland, construction land, and water bodies remained stable. Thus, in the next 20 years, the Xishuangbanna region will undergo further deforestation; simultaneously, given its improved economic development and the rising human demand for resources, the cultivation of cash crops such as rubber and tea will continue to increase, which will add pressure on the region's forests.

**Figure 7.** Land use forecast for Xishuangbanna region for the next 20 years (unit: ha or 104 m2).

#### *3.4. Reference Levels in Xishuangbanna*

According to IPCC's Good Practice Guidance, the source/or sink estimates were determined by multiplying the activity data by a carbon stock coefficient (i.e., emission factor) at two points in time. In this study, the combination of the IPCC method and the land use change model showed that the carbon emissions from the study region obviously increased year by year over the 20 years of this study (Figure 8); the simulated growth trend provides an estimate of 0.35 million t CO2e of annual carbon emissions on average. Simultaneously, the large increase in rubber plantations facilitated a rise in carbon absorption, resulting in average annual carbon sequestration of 0.13 million t CO2e. Although the total amount of carbon sequestration attributable to cultivated land, grassland, shrubland, and tea gardens changed, the overall increase was not large. In general, the total carbon emissions in Xishuangbanna rose year by year during the past two decades. The average annual carbon emissions in the past two decades were estimated to be 0.23 million t CO2e, while the total carbon emissions in the same time period amounted to 4.6866 million t CO2e, indicating an obvious increase.

**Figure 8.** Reference levels in the Xishuangbanna region.

#### **4. Conclusions**

A careful assessment of RLs for REDD+ in Xishuangbanna, China provides significant insights to REDD+ project. The implications that emerge from this study are as follows.

1. We developed a methodological framework to estimate carbon emissions for the REDD+ program in the tropical forests of Xishuangbanna, China. By coupling IPCC's GPG and land use scenario modeling, we could successfully estimate the RLs. Within the framework, the Enhanced Thematic Mapper Plus/Thematic Mapper(ETM+/TM) remote sensing images in the study area were used to interpret land use changes in 1990, 2000, and 2010. The Land Use Dynamic Index was used for the transfer of land use types during the study period to identify that rubber plantations were the main contributor to forest loss in this region. The Markov chain model was used to predict the number of future land use changes and the Multilayer Perceptron model was applied to calculate the distribution location of these changes.

2. According to Paragraph 71 of Decision 1/CP.16, forests RLs are one of the elements to implement REDD+ activities for developing country parties. Moreover, the COP recognizes the importance and necessity of adequate and predictable financial and technology support for developing such RLs. Identifying these RLs is, therefore, a critical step in the provision of financial incentives and/or creation of carbon markets. Furthermore, they guide the design of the REDD+ strategy. In this study for the business-as-usual scenario of the RLs, Xishuangbanna will lose 158,535 ha (158,535 × 104 m2) of forest area in next 20 years, resulting in approximately 0.23 million t (0.23 × 109 kg) CO2e emissions per year. This is due to the improved economic development and the rising human demand for resources, such as the cultivation of rubber and tea.

#### **5. Future Scope**

Estimating carbon emissions based on RLs is a multidisciplinary task. It requires expertise in forestry science, ecological modeling, statistics, remote sensing, and field techniques. Undertaking this exercise is demanding given global geographical diversity, and, thus, building technical capacity to this end is essential. Modeling future emissions based on historical trend rates and understanding the relationships between deforestation patterns and the drivers of deforestation are essential for RL estimation [68–70].

Remote sensing technology using optical sensors is capable of measuring the carbon content of different forest types when supported by field information from, for example, sample plots used to calibrate the technology. Using this methodology, a multitemporal set of remotely sensed data can be used to detect forest changes over time [71–73]. Thus, freely available Landsat images can provide reliable measurements of forest change, especially when complemented with high-resolution satellite imagery from sensors such as QuickBird, which provide data for image analysis training and validation.

**Author Contributions:** Conceptualization, H.L.; Data curation, G.L., C.Z. and K.H.; Formal analysis, K.H., Y.F. and M.X.; Investigation, Y.F., M.X. and R.G.; Methodology, G.L. and C.Z.; Project administration, H.L.; Resources, G.L.; Writing—original draft, H.L. and G.L.; Writing—review and editing, H.L. All authors have reviewed and agreed to the published version of the manuscript.

**Funding:** This study is under the auspices of NSFC42071267, NSFC41371525, Program for Innovative Research Team (in Science and Technology) in University of Henan Province (21IRT-STHN008), Dabieshan National Observation and Research Field Station of Forest Ecosystem at Henan andSYL20060111.

**Acknowledgments:** We are grateful to Quntao Yang's contribution to the research. We are also indebted to Liangyun Liu and Bowen Song, our colleagues at the Aerospace Information Research Institute (AIR) under the Chinese Academy of Sciences (CAS), for their kind help on China's GlobeLand30.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** The accuracies in 1990.

#### **Table A2.** The accuracies in 2000.

**Table A3.** The accuracies in 2010.

#### **References**

1. FAO. *Global Forest Resource Assessment 2005*; FAO: Rome, Italy, 2006.

## *Article* **Assessment of Land Degradation in Semiarid Tanzania—Using Multiscale Remote Sensing Datasets to Support Sustainable Development Goal 15.3**

### **Jonathan Reith 1,2,3,\*, Gohar Ghazaryan 2,4,5, Francis Muthoni <sup>3</sup> and Olena Dubovyk 2,4**

**Abstract:** Monitoring land degradation (LD) to improve the measurement of the sustainable development goal (SDG) 15.3.1 indicator ("proportion of land that is degraded over a total land area") is key to ensure a more sustainable future. Current frameworks rely on default medium-resolution remote sensing datasets available to assess LD and cannot identify subtle changes at the sub-national scale. This study is the first to adapt local datasets in interplay with high-resolution imagery to monitor the extent of LD in the semiarid Kiteto and Kongwa (KK) districts of Tanzania from 2000–2019. It incorporates freely available datasets such as Landsat time series and customized land cover and uses open-source software and cloud-computing. Further, we compared our results of the LD assessment based on the adopted high-resolution data and methodology (AM) with the default medium-resolution data and methodology (DM) suggested by the United Nations Convention to Combat Desertification. According to AM, 16% of the area in KK districts was degraded during 2000–2015, whereas DM revealed total LD on 70% of the area. Furthermore, based on the AM, overall, 27% of the land was degraded from 2000–2019. To achieve LD neutrality until 2030, spatial planning should focus on hotspot areas and implement sustainable land management practices based on these fine resolution results.

**Keywords:** land degradation neutrality; SDG; land productivity; land cover; NDVI; Landsat; vegetation-precipitation relationship; soil organic carbon; Google Earth Engine

#### **1. Introduction**

Land degradation (LD) is defined as the "continuous reduction or loss of the productivity of the land due to a combination of natural and anthropogenic causes" [1]. It is a global problem and affects people, their livelihoods and nature. Studies suggest that up to 3.2 billion people live and depend on degraded lands [2] and that approximately a quarter of the world's lands are affected by LD [3,4]. Poor people, who often rely on agriculture, are most vulnerable to LD [5,6]. Lost ecosystem services due to land use and land cover (LULC) change and LD account for up to USD 10.5 trillion loss per year, which is about a sixth of the world's gross domestic product (GDP) [7]. Furthermore, biodiversity is declining globally, with tremendous losses in sub-Saharan Africa because of LD [6]. Projections suggest that lower productivity in the face of climate change will drive LULC change globally. Moreover, the population growth, combined with a changing diet, will have an enormous influence on agriculture and thus LD [8]. It is for these reasons that the world community introduced the sustainable development goal (SDG) 15.3, which aims to

**Citation:** Reith, J.; Ghazaryan, G.; Muthoni, F.; Dubovyk, O. Assessment of Land Degradation in Semiarid Tanzania—Using Multiscale Remote Sensing Datasets to Support Sustainable Development Goal 15.3. *Remote Sens.* **2021**, *13*, 1754. https:// doi.org/10.3390/rs13091754

Academic Editor: Elias Symeonakis

Received: 31 March 2021 Accepted: 26 April 2021 Published: 30 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

"restore degraded land and strive to achieve an LD-neutral world", highlighting the global importance of this issue [9,10].

Tanzania is a hot spot of LD, with more than half its area showing signs of degradation [2,11]. It has the highest annual forest area net loss in East Africa and the fifthhighest worldwide [12]. The cost of LD has been summed up to USD 2.3 billion annually in the first decade of the new millennium [13]. Seventy-five percent of the total labor force, mostly rural people, work and depend on the agricultural sector, which is accountable for about 30% of the GDP [14]. Although the cultivated area increased in the last years, the output per hectare (ha) decreased, both in annual and perennial crops, even though fertilizer consumption quadrupled at the same time [15]. The number of undernourished people is growing and is currently more than 30% [16]. The population is increasing while agricultural productivity is stagnating, and the economic dependency on natural goods is still high. The consequences of this dilemma area persisting pressure on land and, thus, a probable conversion of natural into cultivated land in the coming years. The poor people's food security is also at risk, and in the coming years, in the face of climate change, new insecurities are likely to arise [17]. This holds especially true for the rural semiarid central districts of Kiteto and Kongwa (KK).

Agricultural intensification and sustainable land management (SLM) are keys to halt and reverse LD [18–20]. One major constraint that prevents action is the lack of spatial information on the extent and magnitude of LD [18]. In contrast to the laborious fieldwork, remote sensing offers the unique opportunity to consistently assess vast areas over a long period [2–4]. Unfortunately, the existing LD maps have a coarse spatial resolution and provide inconsistent estimates of the affected area [8]. For example, previous estimates of the extent of LD in Tanzania range from 41% to half of the country [2,3,11]. These variations emanate from differences in definitions of LD, monitoring methods and lack of appropriate data [6,21]. In the course of SDG 15.3 implementation, standard methods for assessing LD were introduced, making reports more comparable.

This new standard methodology, recommended by the United Nations Convention to Combat Desertification (UNCCD), includes the usage of three sub-indicators for the complimentary assessment of LD [22]. The first sub-indicator, land cover (LC), reports changes in vegetation cover. The second, land productivity (LP), captures changes in ecosystem functions. The last, soil organic carbon (SOC), indicates slower changes resulting from biomass alterations [20]. The three sub-indicators are aggregated to form the land degradation indicator. Improvements in one indicator cannot compensate losses in others, as they are complementary and not additive. Thus, the "one-out, all-out" approach is applied whereby even if one indicator shows signs of decline and the others are positive, the land is deemed to be degraded [23].

The recent Tanzanian national LD-neutrality (LDN) report follows these guidelines [24]. However, it only assesses LD for the first ten years of the 21st century and mainly uses global default data with a coarse spatial resolution. The 1 km coarseresolution is inadequate to monitor LD in small mountainous and highly fragmented landscapes, as it may miss out on smaller than pixel size LD areas [25].

Overall, only a few studies have been published on the subject of SDG 15.3.1 monitoring and assessment. Gichenje and Godinho [26], for example, conducted a baseline assessment of the SDG indicator 15.3.1 for the years 1992 to 2015 using the Advanced Very High Resolution Radiometer (8 km, AVHRR) Normalized Difference Vegetation Index (NDVI) time series and the European Space Agency (ESA) Climate Change Initiative (CCI) LC map in Kenya. In Mozambique, Frederique et al. [27] analyzed only the LD sub-indicator LP trend using the Moderate Resolution Imaging Spectroradiometer (250 m, MODIS) NDVI from 2001 to 2016.

However, these studies share the common disadvantage of applying only default methodology and global datasets for national and subnational LD assessments. Though Akinyemi et al. (2020) used a customized 30m resolution LC map to assess the LC subindicator of SDG 15.3.1 in Botswana, this study relied on AVHRR time-series assessment

for the LP sub-indicator. Furthermore, no studies exist in Africa that used high-spatialresolution datasets the assessment of more than one sub-indicator of SDG 15.3.1. Therefore, it is vital to overcome the existing research gaps and use high-resolution spatial data to provide improved information on the SDG 15.3 [28].

In this light, the main aim of our study was to assess the SDG 15.3.1 indicator based on the newly adopted approach based on the higher resolution (compared to default UNCCD datasets) 30m Landsat time series and 30m LC maps and compare our results to the estimates of the SDG 15.3.1 based on the default UNCCD data and methods.

Our study addressed the following research questions:

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study site is situated in Kiteto and Kongwa districts, located in Dodoma and Manyara regions of Central Tanzania, respectively (Figure 1). The elevation ranges between 850 and 2100 m above sea level. The study area has a hot arid steppe climate [29]. The average monthly temperature stays between 19 and 25 ◦C all year, and the precipitation is roughly 600 mm a year, with interannual differences of 500 to 800 mm. Large parts of northern Kiteto and more minor areas of the mountainous region in Kongwa are protected areas for nature and landscape conservation.

#### *2.2. Materials*

The SDG 15.3.1 indicator and its three LDN sub-indicators were computed using the recommended default method (DM) with Trends.Earth [30] and the adapted methods (AM) using high-resolution Landsat (and other) datasets (Table 1).

**Figure 1.** Location of the study area in Central Tanzania (**A**,**B**) and protected areas (**C**).

The DM (LC) map provided by the UNCCD is based on the 300 m ESA CCI LC map (Table 1). The AM utilized 30 m LC maps for 2000–2018 in the study area that the Regional Centre for Mapping of Resource for Development (RCMRD) developed. Both datasets were disaggregated into the six LC classes as defined by Intergovernmental Panel on Climate Change (IPCC), i.e., forestland, grassland, cropland, wetland, urban, and otherlands [31].

The recommended global default dataset uses the MOD-13Q1-coll6 (250 m) MODIS-NDVI products [30]. In contrast, the AM was calculated based on a 30 m resolution NDVI from a combination of Landsat 5, 7 and 8 (Table 1). The Landsat time series were accessed and analyzed using Google Earth Engine [32], based on atmospherically corrected surface reflectance collections (Table 1). The Landsat 5 and 7 data were spectrally harmonized with Landsat 8 series using linear transformation [33]. As a further step to improve the image quality, the fmask was adopted to mask out clouds and cloud shadows [34,35]. Generally, the images with cloud cover scores higher than 80% were removed. Finally, the NDVI was calculated for each image, and then the images of the same admission time were merged and clipped to the extent of the study area. As it is recommended to constrain the observation period to the growing season to reduce the number of irrelevant assets for the computation and enhance the quality of the time series [22], we used the imagery from November to June. When using Trends.Earth, there is no possibility to apply the computation to the growing season, so the DM uses the whole calendar year. In order to integrate the rainfall information, data from Climate Hazards Group InfraRed Precipitation with Station (CHIRPS) were used (Table 1).

The SOC metrics were derived from the SoilGrids250m dataset [41] for the DM and the AM, as there is no national SOC database for Tanzania. SOC is measured at a depth of 30 cm and is stated as mass per area (e.g., tons per hectare (t/ha)) [22].

#### *2.3. Methods*

The calculation of the SDG 15.3.1 indicator is based on the "one out, all-out" approach (Ref. [23] and Figure 2). The three LD sub-indicators (LC change, LP decline and loss of SOC) are estimated, and if one indicator signals degradation, the LD indicator will reflect this as well. A baseline is needed to compare the progress of LDN. The baseline year (t0) was set to be 2015 and is computed as the average of the period leading up to t0 (2000–2015). The indicators are then remeasured in regular time intervals leading to 2030, and change is used to monitor the progress to accomplish LDN [20].

**Figure 2.** Steps to derive the sustainable development goal (SDG) indicator 15.3.1 from the sub-indicators. I represents Improvement, S represents Stable and D represents degraded (based on [31]).

To calculate the indicator for the reporting year 2019 (t1), it is necessary first to assess the baseline util t0 and then calculate the change from the baseline to t1 (Figure 2). As a final step, combine both results. The details of the calculation of each indicator are explained in the following section. The three LD sub-indicators were created from satellite images using cloud-based geospatial computing. The indicators were calculated using Trends.Earth [30] and Google Earth Engine [32] for the DM and AM, respectively. As Trends.Earth currently only enables the computation for the baseline period (BP), the DM is only available from 2000 to 2015.

#### 2.3.1. Sub-Indicator 1: Land Cover Transitions and Degradation

The first SDG 15.3.1 indicator is the LC change. To assess the LC degradation, the transitions between 2000–2015 and 2015–2018 were analyzed for the baseline and the first monitoring period (MP), respectively. To determine whether changes from one LC class to another are interpreted as degradation, a change matrix can help visualize the transitions (Table 2) based on the Good Practice Guidance by the UNCCD [31]. It is recommended to adopt this matrix for the national context. Therefore, transitions from grasslands to croplands were not considered LD for the AM to avoid tradeoff between ecosystems and food security and between nomadic and sedentary living.

#### 2.3.2. Sub-Indicator 2: Loss of Land Productivity

LP is described as "the biological productive capacity of the land". It is closely associated with net primary productivity [42], which can be measured directly with earth observation methods [22]. NDVI is a widely used index detecting LP [26,43,44]. The LP sub-indicator consists of three distinct components, namely trend, state and performance.

The LP trend component measures the trajectory of change in productivity over time. It is calculated at the pixel level using linear regression and the Mann Kendall significance test [22,45,46]. Positive and negative changes in NDVI indicate increasing and decreasing productivity associated with vegetation recovery and degradation, respectively. The eight most recent years of data were used to create a new distinct and significant time series that is more responsive to present land conditions. Further, following [47], we accounted

for the effect of rainfall variability on vegetation productivity trends by using the rain use efficiency (RUE) method.

The LP state component represents recent changes in LP compared to the BP. The yearly NDVI mean images of the shortened BP (2000–2012) were normalized and assigned to classes from 1 to 10 based on their percentiles. To avoid annual fluctuations, contemporary values of the three-year anteceding t0 and t1 were classified in this scheme. Areas with a reduction of two or more classes were classified as degraded, while the rise by two categories was interpreted as an improvement [31].

The LP performance component examines local productivity compared to similar ecoregions defined by the unique combination of SoilGrids [41], soil taxonomy great groups and LC classes (Table 1). The 90th percentile in each ecoregion was calculated as a proxy for the maximum productivity level. The LP performance was then calculated based on the ratio of the observed mean NDVI value per pixel and the NDVImax (90th). Values below 0.5 indicate regions where the LP is low and LD may prevail [31].

The overall LP sub-indicator is calculated based on the three components mentioned earlier. As the LP trend is based on a statistically significant test, it is most influential, and its status determines LP degradation. Only if both LP status and LP performance show negative results, does the LP indicator also show degradation [22]. If only the LP state component shows degradation, this could indicate "early signs of decline" because the other indicators may not have detected the most recent LD. Further, if only performance shows degradation, there is no temporal trend, and the land is classified as "stable but stressed" [22]. In contrast to the Good Practice Guidance by UNCCD, Trends.Earth (DM) also incorporates the "early signs of decline" state component into the LP degradation [30].

**Table 2.** Land cover transition matrix (2000–2015) based on the adapted methods (AM). Green, beige and brown colors indicate improving, stable and declining conditions of land cover categories, respectively. The area in km2 and the possible cause of the land cover transition are indicated in the matrix. The change is based on the high-resolution land cover dataset.

#### 2.3.3. Sub-Indicator 3: Degradation of Soil Organic Carbon

The Good Practice Guidance for the SOC sub-indicator is based on the maximum equilibrium SOC content at a location that is controlled by environmental factors such as rainfall, evaporation, solar radiation, and temperature [22]. The content can change based on three distinct change factors: First, the land-use factor represents SOC stock changes based on the type of land use. Second, the management factor reflects the management practice of the land use (e.g., grazing intensity on grasslands). Third, the input factor represents the different amounts of carbon input into the soil [22,48,49]. While the LULC change factor can be used with LC as a proxy, there are presently no sufficient datasets available to provide information about the management or the input for the other two indicators. Thus, the only indicator to assess SOC changes is the second LD indicator LC change [22].

#### **3. Results**

Three sub-indicators, namely LC transitions, LP decline, and SOC loss, were estimated to derive the SDG 15.3.1 indicator using the default and adapted methods. The patterns of each sub-indicator based on DM and AM are described in the following sections starting with the BP from 200 to 2015 for both DM and AM. The first monitoring period from 2015 to 2019 is only assessed using the AM, as the data necessary for this period are currently not available in Trends.Earth.

#### *3.1. Sub-Indicator 1: Land Cover Transitions and Degradation*

According to the DM based on the medium-resolution 300 m LC maps, over 99% of the study area remained stable in the BP (2000–2015) (Table A1). Urban areas covering less than 0.1% of the study area experienced the highest relative expansion (56%). The forestlands were the only other LC class that increased in area significantly (4.4%) in the BP.

In contrast to the DM, the AM with high-resolution (30 m) LC data revealed that 6.7% of the total area changed to a less desirable LC class, signifying LD, and only 2.3% of analyzed areas improved. The area of (semi)natural LC, such as forestlands (−19%), grasslands (−6.6%) and wetlands (0.1%), mostly declined, whereas the croplands recorded the highest spatial gain (24.2%) (Figure 3 and Table 2).

**Figure 3.** Sankey plot describing the land cover transitions between the years 2000, 2015 and 2018 using high-resolution land cover data. Bands represent the actual proportion of land that changed class over time.

The trend observed in the BP continued in the first years of the MP (Figure 3 and Table A2). Overall, from 2015 to 2018, 3.3% of the total area was degraded during the MP, while 1.2% of the area changed to a more desirable LC. Grass- and forestlands continued to decline by 3 to 9%, respectively, while anthropogenic(-influenced) covers such as cropland and urban areas expanded. Compared to about 3000 ha forests lost per year (a) in the BP, the rate doubled to 6000 ha/a in the MP. Similarly, the changes in croplands increased from 6000 ha/a in 2000–2015 to 7500 ha/a in 2015–2018.

#### *3.2. Sub-Indicator 2: Loss of Land Productivity*