Combining Earth Observations, Cloud Computing, and Expert Knowledge to Inform National Level Degradation Assessments in Support of the 2030 Development Agenda

Abstract

Monitoring progress towards the 2030 Development Agenda requires the combination of traditional and new data sources in innovative workflows to maximize the generation of relevant information. We present the results of a participatory and data-driven land degradation assessment process at a national scale, which includes use of earth observation (EO) data, cloud computing, and expert knowledge for Argentina. Six different primary productivity trend maps were produced from a time series of the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) dataset (2000–2018), including the most widely used trajectory approach and five alternative methods, which include information on the timing and magnitude of the changes. To identify the land productivity trend map which best represented ground conditions, an online application was developed, allowing 190 experts to choose the most representative result for their region of expertise nationwide. Additionally, the ability to detect decreases in land productivity of each method was assessed in 43,614 plots where deforestation had been recorded. The widely used trajectory indicator was the one selected by most experts as better reflecting changes in land condition. When comparing indicators’ performance to identify deforestation-driven reductions in productivity, the Step-Wise Approach Trend Index (SWATI), which integrates short- and long-term trends, was the one which performed the best. On average, decreases of land productivity indicate that 20% of the Argentine territory has experienced degradation processes between 2000 and 2018. The participatory data generation and verification workflow developed and tested here represents an innovative low cost, simple, and fast way to validate maps of vegetation trends and other EO-derived indicators, supporting the monitoring of progress towards land degradation neutrality by 2030.

1. Introduction

Land degradation is the consistent loss of ecosystem functionality due to human and natural processes [1]. It affects biodiversity and ecosystem services, exacerbates climate change, and ultimately impacts the well-being and livelihoods of 1.5 billion people globally [1,2]. In response to the dramatic socio-environmental challenges humanity is facing, the United Nations (UN) adopted the 2030 Agenda for Sustainable Development. A total of 17 Sustainable Development Goals (SDGs) came into effect on January 2016 and are expected to guide social, economic, and environmental policy and investment over the next 15 years. In particular, SDG 15 promotes “Life on Land” and SDG Target 15.3 focuses on achieving land degradation neutrality (LDN) by 2030, by combating desertification, restoring degraded land and soil, including land affected by desertification, drought, and floods, and striving to achieve a land degradation neutral world. LDN is defined as “a state whereby the amount and quality of land resources necessary to support ecosystem functions and services and enhance food security remain stable or increase within specified temporal and spatial scales and ecosystems” [3].

The achievement of LDN is a complex process that includes numerous political, technical, and financial aspects. The United Nations Convention to Combat Desertification (UNCCD) LDN Technical Guide proposes a step-wise approach to define LDN targets and identify measures to achieve them [4], for which the LDN scientific conceptual framework flags the importance of both the assessment of a baseline and the monitoring of indicators [3]. A set of indicators are proposed by the LDN framework for country parties to inform on the proportion of land that is degraded over total land area (SDG Indicator 15.3.1) and to report on progress towards the UNCCD objectives. These three indicators are: trends in land cover (SO1-1), trends in land productivity or functioning of the land (SO1-2), and trends in carbon stocks above and below ground (SO1-3). Reporting on these land-based indicators is encouraged by being supported by comparable and standardized nationally validated data sources.

National data on the three indicators can be collected through existing sources (e.g., databases, maps, reports), including participatory inventories on land management systems, as well as remote sensing data collected at the national level. Regional and global datasets derived from global earth observation (EO) systems play an important role in the absence of, to complement, or to enhance national official data sources, particularly for assessing changes in land productivity and land cover dynamics [5]. Since just three indicators cannot fully capture the complexity of land degradation (i.e., its degree and drivers), countries are strongly encouraged to use other relevant national or sub-national indicators, data, and information to strengthen their interpretation, as well as participatory processes to validate results based on EO data. Integration of EO data into participatory processes that include local knowledge is crucial in the process of LDN [3,4,5,6].

Particularly, the trends in the land productivity indicator (SO1-2) refer to changes in the total above-ground net primary production (NPP), defined as the energy fixed by plants minus their respiration, which translates into the rate of biomass accumulation that delivers a suite of ecosystem services. This indicator points to changes in the health and productive capacity of the land and reflects the net effects of changes in ecosystem functioning due to changes in plant and biomass growth, where declining trends are often a defining characteristic of land degradation [7]. Understanding changes in the productive capacity of the land is critical for assessing the impact of land management interventions, its long-term sustainability, and the climate-derived impacts which could affect ecosystem resilience and human livelihoods. NPP is usually estimated by means of EO from known correlations between the fraction of absorbed photosynthetically active radiation (fAPAR) and plant growth vigor and biomass [8]. Various vegetation indices can be calculated from satellite data to be used as a proxy for land productivity [9], with the Normalized Difference Vegetation Index (NDVI) being the most frequently used [10,11,12]. There are also currently numerous methods that can be used to analyze and extract insights from time series image datasets [13,14,15,16,17,18,19]. For example, the UNCCD Good Practice Guidance (GPG) [9] suggests three metrics derived from the remotely sensed data to estimate changes in land productivity (trend, state, and performance), and tools, such as Trends.Earth, have been developed to support the calculation of such metrics globally [20]. Additionally, many annual metrics linked to key phenological events (i.e., onset of greenness, time of peak NDVI, maximum NDVI, rate of green-up, rate of senescence and integrated NDVI) can be derived from the NDVI time series. These metrics can be later used to calculate annual trends and can also be combined to create new indices. The choice of the most representative metric to asses land productivity dynamics represents a challenge, especially in large countries, with contrasting environmental conditions, ecosystems, and land uses. Different calculations based on the same EO data can produce different results [21], highlighting the importance of integrating EO data with other sources of information, such as experts’ knowledge through participative processes [6,22]. Accordingly, the UNCCD LDN reporting framework encourages countries to assess direct and/or indirect drivers of land degradation through participatory processes adapted to the national and local context.

In this paper, we present the results of a participatory and data-driven land degradation baseline assessment process at a national scale, which includes use of EO data, cloud computing, and expert knowledge for Argentina. Our main goals included (1) obtaining primary productivity trend maps at a national level from the time series of the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI data (2000–2018), including six alternative methods to estimate this indicator; (2) developing and testing an online application that allows experts to remotely choose the most representative result in their region of expertise and identify the main drivers of change; and (3) comparing the ability to detect decreases in primary productivity of each indicator in areas of known deforestation.

2. Materials and Methods

2.1. Study Area

With a continental area of 2,780,400 km², Argentina is the eighth-largest country in the world. Drylands represent 55% of the country’s land area, and, to varying degrees, all of them present degradation problems [23]. Due to the diversity of environmental conditions, ecosystem types span from subtropical rainforest in the north, to the Patagonia desert sub-Antarctic areas in the south, with regions of humid pampas in the center and the east, as well as a vast semiarid to hyper-arid diagonal region towards the Andes. Land degradation is a major challenge in the country and the magnitude of the economic losses caused by these processes represent about 16% of the national GDP [24].

Argentina ratified the UNCCD in 1996 and the Secretariat of Environment and Sustainable Development prepared the National Action Plan to coordinate all major sustainable land management (SLM) projects and programs to implement the UNCCD objectives. Argentina was one of the six pilot countries selected by the Food and Agriculture Organization of the United Nations (FAO) for studying Land Degradation Assessment in Drylands (LADA). As a follow-up to LADA, the National Observatory of Land Degradation and Desertification was established to help formulate policies and strategies for controlling and mitigating land degradation and desertification. The tools and processes presented in this work were developed in the context of the UNCCD fifth reporting cycle (PRAIS 3) and the LDN Target Setting Programme of the Global Mechanism.

2.2. Time Series of Earth Observation Data

To calculate the different primary productivity trend indicators for the period 2000–2018, a time series of the NDVI dataset from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1) Version 6 was used. These data are generated by the United States National Aeronautics and Space Administration’s Land Processes Distributed Active Archive Center (NASA LP DAAC) every 16 days at 250 m spatial resolution, producing 23 composites per year [25]. The MODIS NDVI product is computed from atmospherically corrected bi-directional surface reflectance that has been masked for water, clouds, heavy aerosols, and cloud shadows. The algorithm chooses the best available pixel value from all the acquisitions from the 16 day period. To capture the full annual cycle of vegetation in the southern hemisphere, the phenological year was considered, which begins on July 1 and ends on June 30 [26].

To summarize MODIS NDVI data from the year 2000/2001 to 2017/2018 (414 composites corresponding to 18 years), two different annual metrics were calculated and, therefore, two time series of annual data from 2000 to 2018 were obtained for every pixel:

Annual Mean (AM): The AM was calculated as the arithmetic mean [27] of the NDVI values within each year, by first adding the 23 NDVI values and then dividing the sum by 23. This mean annual NDVI has the same behavior as the accumulated NDVI normally associated to annual NPP, but with the advantage of keeping values in the NDVI value range;

Annual Ecosystem Services Productivity Index (ESPI): The ESPI is the product between the annual mean and the difference between 1 and the coefficient of variation [27]: ESPI = AM × (1 − CV) [28]. The ESPI does not directly measure ecosystem services but is an index that merges two attributes of the seasonal dynamics of the NDVI; the annual mean and the intra-annual variability (CV). The inclusion of the seasonality helps to better differentiate pasture/crop transitions, which in many cases is difficult due to the fact that both land covers may present the same mean annual NDVI. Intra-annual variability in crops is much higher than in pastures, resulting in a lower ESPI (see Appendix

Table A1. Value and trend intensity categorization for each method and study case.

Indicator	Value	Trend Intensity		Value	Trend Intensity
	CASE 1			CASE 3
LTT-AM	−4%	0	No trend	8%	0	No trend
LTT-ESPI	−58%	3	Strong negative	62%	3	Strong positive
SWATI-AM	−1	1	Light negative	0	0	No trend
SWATI-ESPI	−3	2	Moderate negative	2	2	Moderate positive
ToC-BFAST	−1	1	Light negative	1	1	Light positive
SSWATI-ESPI	−39%	2	Moderate negative	52%	3	Strong positive
	CASE 2			CASE 4
LTT-AM	43%	2	Moderate positive	−11%	0	No trend
LTT-ESPI	64%	3	Strong positive	23%	0	No trend
SWATI-AM	3	2	Moderate positive	1	1	Light positive
SWATI-ESPI	3	2	Moderate positive	1	1	Light positive
ToC-BFAST	0	0	No trend	0	0	No trend
SSWATI-ESPI	60%	3	Strong positive	106%	3	Strong positive

and Figure A1 for single pixel examples illustrating this).

2.3. Trends in Land Productivity

With the objective of comparing different methodologies for characterizing land productivity trends, alternative methods were used:

Long-Term Trend (LTT): This method is based on the trajectory indicator as described by the UNCCD GPG [9]. The trend is obtained by computing a Mann–Kendall non-parametric significance test [29,30] and considering only as significant changes those that show a p-value ≤ 0.05. Additionally, a modification was introduced in which, for all pixels with a significant trend, the Thiel–Sen slope [31] was calculated as a measure of the trend intensity. This intensity (rate of change) was expressed as a percentage of the time series initial value, to compare among different pixels.

Step-Wise Approach Trend Index (SWATI): This method was developed to capture the inter-annual variability present in different sections of the time series. For the calculation, the time series is divided into sections (or steps), inside of which the significance of the trend is evaluated with the Mann–Kendall non-parametric test. Significantly negative trends are assigned a score of -1, a score of 0 is set for non-significant trends, and a score of 1 is given to significant positive trends. This method was designed for detecting non-linear behaviors or short term trends that the LTT cannot capture. For these analyses, the step duration was set to 5 years, in accordance with the planned monitoring cycles, leaving a baseline period of 8 years at the start of the time series: Baseline Step (BS: 2000/2001 to 2007/2008), Step 2 (S2: 2008/2009 to 2012/2013), and Step 3 (S3: 2013/2014 to 2017/2018) and the long-term trend (LTT: 2000/2001 to 2017/2018). The indicator is computed as SWATI = BS + S2 + S3 + 2*LTT. The resulting values vary, in this case between -5 and +5, indicating different levels of intensity trend (negative and positive, respectively).

Slope Step-Wise Approach Trend Index (SSWATI): This method follows the same principle of the SWATI by considering fixed time intervals, but instead of assigning scores to each step, the magnitude of the change inside that step is taken into account. The magnitude of each time step is computed as the slope (Thiel–Sen slope) multiplied by the number of years in the step. This is only calculated when the trend is significant, otherwise the value of change is zero. The calculation using the final magnitudes (m) is performed as SSWATI = mBS + mS2 + mS3.

Type of change using BFAST (ToC-BFAST): This method uses the 16-day MODIS MOD13Q1 product to directly analyze NDVI using the breaks for additive seasonal and trend (BFAST) procedure [32,33]. BFAST is an iterative algorithm that combines the decomposition of the time series into seasonal, trend, and residual components and allows the detection of trend changes within the time series, while accounting for seasonality. This indicator was provided by de Abelleyra et al. [34], who searched the time series not only for the LTT, but also for a single break that produces a two-part trend. When present, the year of the breakpoint was registered, as well as the type of trend before and after the breakpoint. According to this behavior, eight classes of types of change are described in their work. De Abelleyra et al. [35] did not include the high Andean mountain chain that runs along the west border of Argentina from north to south, to discard pixels that are frequently covered by snow.

LTT and SWATI methods were calculated both with the annual mean NDVI time series (LTT-AM and SWATI-AM) and with the ESPI time series (LTT-ESPI and SWATI-ESPI), SSWATI was calculated with the ESPI (SSWATI-ESPI), and the ToC-BFAST was calculated with the original 16-day NDVI product time series. Thus, a total of six primary productivity trend maps at a national level were obtained.

2.4. Categorization of Trend Intensity

For the purposes of reporting on SDG Indicator 15.3.1, it is not necessary to quantify the magnitude of change in productivity in biomass units of net primary productivity (NPP); rather, it is important to determine whether productivity is increasing (greening: positive trend), decreasing (degradation: negative trend), or if it remains stable (no change: no significant trend) for the land unit over time [36]. However, some indication of the magnitude of the change (for example, identifying if the rate of change represents a strong, medium, or light variation) is of much importance in the LDN scheme. The LDN framework is based on hierarchical responses, establishing different actions depending on the level of degradation [3]. To compare the different primary productivity trend maps, results were classified in seven categories of trend intensity (TI) according to the magnitude and type of trend, ranging from −3 (strong negative trend) to 3 (strong positive trend). Negative TIs correspond to negative trends (strong, medium, and light), which indicate a decrease in productivity, whereas positive TIs correspond to positive trends, that indicate that productivity is increasing. For the LTT methodology, the magnitude of change was calculated as the relation between the total variation (final estimated value – initial estimated value) and the initial value, expressed as percentage, where the initial estimated value corresponds to the intercept of the adjusted linear model at year 2000. A TI of zero indicates that no changes in productivity were detected with a significant level. In this way, the TI provides seven qualitative classes of land productivity trends for the 2000−2018 period, as shown in

Table 1. Trend intensity (TI) categorization for the different methods. When changes are expressed as percentages, they are referred to the initial value (intercept at year 2000 of the estimated linear model). LTT: Long-Term Trend, AM: Annual Mean; ESPI: Annual Ecosystem Services Productivity Index; SWATI: Step-Wise Approach Trend Index; ToC-BFAST: Type of change using BFAST (breaks for additive seasonal and trend); SSWATI: Slope Step-Wise Approach Trend Index.

TI	Description	LTT (AM and ESPI)	SWATI (AM and ESPI)	ToC-BFAST	SSWATI-ESPI
−3	Strong Negative trend	Decrease of at least 50%	SWATI = −4 or −5	ToC = 1, 3, or 5 and decrease of at least 50%	SSWATI < 0, decrease of at least 50%
−2	Moderate Negative trend	Decrease between 25% and 50%	SWATI = −2 or −3	ToC = 1, 3, or 5 and decrease between 25% and 50%	SSWATI < 0, between 25% and 50%
−1	Light Negative trend	Decrease of less than 25%	SWATI = −1	ToC = 1, 3, or 5 and decrease of less than 25%	SSWATI < 0, between 25% and 10%
0	No Trend	No significant slope (p > 0.05)	SWATI = 0	ToC = 7 or 8 and increase	Decrease or increase lower than 10%
1	Light Positive trend	Increase of up to 25%	SWATI = 1	ToC = 2, 4, or 6 and increase of up to 25%	SSWATI > 0, between 10% and 25%
2	Moderate Positive trend	Increase between 25% and 50%	SWATI = 2 or 3	ToC = 2, 4, or 6 and increase between 25% and 50%	SSWATI > 0, increase between 25% and 50%
3	Strong Positive trend	Increase of at least 50%	SWATI = 4 or 5	ToC = 2, 4, or 6 and increase of at least 50%	SSWATI > 0, increase of at least 50%

.

2.5. Online Application for Expert Data Collection

Expert knowledge can provide a source of information to help validate and understand the underlying causes of the productivity trend in the territory. In large countries, such as Argentina, gathering expert knowledge for all areas represents a challenge. To gather expert data at a low cost and in a simple and fast way, a web-based tool was developed using open and free software. A JavaScript interface was built to connect the experts via a web page to a GeoServer application [37] and a PostgreSQL database. Each expert was given login credentials (username and password) and a tutorial explaining how to use the online application. After logging in, the interface shows a Google Earth high resolution satellite image of Argentina as a base map and the polygons corresponding to the vegetation units of the country [38]. By clicking on the map, experts could choose the polygon (vegetation unit) to verify, and then they were presented with a panel in which they could visualize each primary productivity trend map, as shown in Figure 1. After viewing the six maps, experts were asked to select the one that best represented the spatial pattern of land productivity trends between 2000 and 2018 in their region of expertise. They were also asked to select, in order of importance, three causes for positive trends and three causes for negative trends. A list of 190 experts in the subject was provided by different institutions, such as The National Council of Science and Technology (CONICET) and the National Institute of Agropecuarian Technology (INTA), which were invited to participate.

2.6. Comparison of Indicators’ Ability to Detect Decreases in Primary Productivity in Plots with Forest Loss

Forest loss is a type of land cover change that generates, in most cases, a decrease in productivity [39,40]. To compare the ability of each method to detect decreases in primary productivity in plots where deforestation had been recorded, we used a dataset produced by Argentina’s Management Unit of the Forest Assessment System (UMSEF) of the National Forest Agency. This dataset includes areas where loss of forest lands and loss of other woodlands occurred (at least 80% of forest cover was lost). Forest lands correspond to land covered by trees higher than 7 m and with a canopy cover of more than 20%. Other woodlands correspond to land with a canopy between 5% and 20% of native trees higher than 7 m or lands with a canopy cover of more than 20% of trees or shrubs shorter than 7 m but higher than 0.5 m, including gallery forests, palm groves, reed beds, and shrublands.

The dataset included 43,614 plots of varying sizes (between 5 and 1000 ha) where deforestation occurred during the studied period, located in four types of forests: subtropical dry forests (Chaco), montane moist forests (Yungas), temperate shrublands (Espinal), and subtropical moist forests (Atlantic Forest). Since forest loss is usually associated with a decrease in primary productivity, we expected negative trends in the deforested plots. The ability of the alternative methods to detect decreases in primary productivity due to deforestation was compared via a linear model where the dependent variable was the trend intensity spatial mean of the deforested plots. In addition to the methods (LTT-AM, LTT-ESPI, SWATI-AM, SWATI-LTT, ToC-BFAST, and SSWATI-ESPI), the type of land (forest land and other woodlands) and the regions (Chaco, Yungas, Espinal, and Atlantic Forest) were included as fixed effects in the model, as well as the interactions among them. When differences were significant, a DGC a posteriori test was performed [41]. We considered that those indicators with lower TIs performed better at detecting decreases in primary productivity due to deforestation.

3. Results

3.1. Trends in Land Productivity

Land degradation has historically been considered a major environmental problem in Argentina, and results found by our analysis confirm that perception. On average, decreases of land productivity indicate that 20% of the Argentine territory has experienced degradation processes between 2000 and 2018. Single model estimates varied between 9% and 35.5% (LTT-AM and ToC-BFAST). Spatially, most models agree with degradation occurring in northern Patagonia and the Dry Chaco, as shown in Figure 2, whereas visible differences between the ToC-BFAST and the rest of the indicators were found, as shown in

Table A2. Percentage of land with each trend intensity (TI) category for each indicator in all ecoregions and the country total area. The average among six indicators is also shown.

Ecoregion	TI	LTT-AM	LTT-ESPI	SWATI-AM	SWATI-ESPI	ToC-BFAST	SSWATI-ESPI	Average
Argentina	−3	0.09	1.68	1.04	0.93	8.97	3.05	2.6
	−2	2.08	4.45	8.57	9.7	8.96	5.69	6.6
	−1	6.73	4.47	13.85	14.11	17.42	6.58	10.5
	0	71.84	73.21	33.17	36.17	40.77	58.79	52.3
	1	12.62	7.72	21.59	20.75	4.22	9.75	12.8
	2	3.57	3.89	17.38	14.85	8.27	9.15	9.5
	3	1.84	3.35	3.18	2.28	4.45	5.78	3.5
	ND	1.21	1.21	1.21	1.21	6.94	1.21	2.2
Altoandina (Montane Grasslands and Shrublands)	−3	1.8	1	6.5	8.9	4.3	2.5	3.8
	−2	0.2	0.3	2.2	1.9	2.7	1.7	1.2
	−1	0.1	0.1	0.1	0.1	9.6	2.8	1.7
	0	53	58.3	31.1	34.7	19.7	67.3	32.8
	1	17.7	9.1	14.5	14.4	2.1	5.8	9.6
	2	8.5	6.2	38.9	35.7	3.4	5.5	15.5
	3	16.8	23.1	4.9	2.5	6.7	12.4	9.0
	ND	1.9	1.9	1.9	1.9	51.5	1.9	9.9
Chaco (Subtropical Grasslands, Savannas, and Shrublands)	−3	8.5	3.7	18.9	14.2	25.2	8.3	11.8
	−2	4.3	5.6	12.4	10.9	9.3	6.7	7.1
	−1	0.2	3.2	2	1.2	5.8	3.7	2.1
	0	70.8	72.7	32	32.6	37.2	54.4	40.9
	1	12.7	9.2	16.9	23.2	5.1	11.8	11.2
	2	2.1	3.7	14.8	14.8	14.2	10.9	8.3
	3	0.6	1.2	2.3	2.3	0.6	3.5	1.2
	ND	0.8	0.8	0.8	0.8	2.4	0.8	0.9
Yungas (Subtropical Montane Moist Broadleaf Forests)	−3	6.9	3.6	13.6	13.5	35.8	6.2	12.2
	−2	1.8	2.5	9.4	6.4	10.4	4.7	5.1
	−1	0.1	1.1	1.3	0.6	3.5	1.8	1.1
	0	75.9	82	33.9	42.5	37.1	63.9	45.2
	1	14.4	8.9	23.9	23.7	4.7	14.7	12.6
	2	0.6	1.3	15.5	11.6	7	7.1	6.0
	3	0.2	0.6	2.3	1.7	0.6	1.5	0.9
	ND	0.1	0.1	0.1	0.1	1	0.1	0.2
Espinal (Temperate Grasslands, Savannas, and Shrublands)	−3	12	3.8	15.6	14	23.2	6.6	12.5
	−2	3	10.7	13.6	14.6	11.9	6.9	10.1
	−1	0.1	3.1	1.4	0.9	8.3	4.5	3.1
	0	74.2	73	34.7	31.4	41.5	50.8	50.9
	1	5.8	1.9	20.5	24.4	4.7	8	10.9
	2	2.5	3.5	9.2	9.5	7.7	11.6	7.3
	3	0.4	1.9	3	3.1	0.3	9.7	3.1
	ND	2.1	2.1	2.1	2.1	2.4	2.1	2.2
Monte (Temperate Grasslands, Savannas, and Shrublands)	−3	1.3	1.8	15.6	16.7	19.2	6.3	10.2
	−2	0.5	0.9	3.6	4	8.9	6	4.0
	−1	0	0.1	0.5	0.6	13.5	1.7	2.7
	0	76.9	77.5	33.1	35.7	40.6	57.9	53.6
	1	10.9	7.9	23.3	21.8	2.2	7.7	12.3
	2	7.3	6	19.1	17.2	6.7	10.3	11.1
	3	2.3	5.2	4	3.2	6.2	9.4	5.1
	ND	0.7	0.7	0.7	0.7	2.8	0.7	1.1
Pampas (Temperate Grasslands, Savannas, and Shrublands)	−3	12.2	4.8	12.4	13.8	17.8	6	11.2
	−2	3.6	11	13.5	16.7	12.4	7	10.7
	−1	0.1	4.1	0.9	1.2	6.8	5.5	3.1
	0	73.1	70.8	34.8	32.3	46.1	52.8	51.7
	1	6.8	2.8	23.5	22.8	5.7	8.8	11.7
	2	2.1	3.2	11.7	10.4	9.5	11.4	8.1
	3	0.9	2.2	2.1	1.9	0.5	7.5	2.5
	ND	1	1	1	1	1.2	1	1.0
Paranense (Subtropical Moist Broadleaf Forests)	−3	6.8	5.4	14.9	14.6	17.7	5.5	10.8
	−2	1	3.5	7.3	9.5	5.1	6.2	5.4
	−1	0.5	1.6	1.1	1	7	3	2.4
	0	55.3	60.4	26.2	29.9	37.7	61.8	45.2
	1	27	18.8	14.1	15.5	7.5	7.8	15.1
	2	1.6	2	26.2	20.7	19.4	5.4	12.6
	3	1	1.5	3.5	2	0.9	3.6	2.1
	ND	6.8	6.8	6.8	6.8	4.6	6.8	6.4
Patagonia (Temperate Grasslands, Savannas, and Shrublands)	−3	5.5	8.4	10.5	12.9	8.3	6.6	8.7
	−2	0.7	1.2	6.1	9.5	7.2	4.1	4.8
	−1	0	0	0.7	0.8	8	1.4	1.8
	0	75.7	77.8	34.6	45.5	49.8	67.1	58.4
	1	13.4	6.2	26.9	18	4.1	10.8	13.2
	2	2.5	3.1	16.7	10.7	4.6	6.3	7.3
	3	1	1.9	3.3	1.2	9.8	2.4	3.3
	ND	1.3	1.3	1.3	1.3	8.2	1.3	2.5
Ecotone Monte- Patagonia (Temperate Grasslands, Savannas, and Shrublands)	−3	10.2	17	31.1	31.2	17.1	14.7	20.2
	−2	2	4.5	17.8	20.3	13.3	21.2	13.2
	−1	0.1	0.1	6.2	7.9	29.4	8.3	8.7
	0	86	77.2	29.1	28.3	31.2	44.6	49.4
	1	0.6	0.2	12.8	10.1	0.6	4.6	4.8
	2	0.3	0.1	2.1	1.5	2.5	5.1	1.9
	3	0.1	0.1	0.3	0.1	4.9	0.9	1.1
	ND	0.7	0.7	0.7	0.7	1	0.7	0.8
Prepuna (Montane Grasslands and Shrublands)	−3	1.3	1.4	6.3	17.1	24.2	9.5	10.0
	−2	0.2	0.2	1.5	1.9	8.6	2.9	2.6
	−1	0	0	0.2	0.2	19.8	0.3	3.4
	0	60.2	62.9	30.5	31.5	24.2	68.4	46.3
	1	30.8	30.2	23.4	17.7	3.2	12.7	19.7
	2	6.9	4.9	33.2	28.6	2.7	5.4	13.6
	3	0.5	0.4	4.9	3.1	7.4	0.8	2.9
	ND	0	0	0	0	9.8	0	1.6
Puna (Montane Grasslands and Shrublands)	−3	1.1	1.2	6.1	12.6	14.5	4.6	6.7
	−2	0.1	0.2	1.5	2.2	6.1	1.9	2.0
	−1	0.1	0.1	0.1	0.1	18.9	2.3	3.6
	0	52	63.1	28	36	28.1	69.4	46.1
	1	37.9	28.3	17.8	14.5	4.3	14.4	19.5
	2	6.1	4	39.2	29.8	3	5	14.5
	3	2.5	3	7.3	4.6	14.4	2.2	5.7
	ND	0.1	0.1	0.1	0.1	10.8	0.1	1.9
SubAntartic (Temperate Grasslands, Savannas, and Shrublands)	−3	3.9	1.7	8	9.4	6.5	3.8	5.6
	−2	0.4	0.6	3.8	2.5	3.8	2.7	2.3
	−1	0.1	0.1	0.1	0.1	0.9	1.1	0.4
	0	70	73.4	36.9	41.1	19	62.2	50.4
	1	15.6	10.5	22.3	20.7	1.9	8.1	13.2
	2	2.5	3.8	20	18.7	6.7	8	10.0
	3	2.1	4.5	3.5	2.1	0.2	8.7	3.5
	ND	5.4	5.4	5.4	5.4	61.1	5.4	14.7

. For all methods, the most frequent category of trend intensity (TI) was “No Trend”, however there was a large range in the percentages categorized as such (between 33% and 73%), as shown in Figure 3. The percentage of land with positive trends was 26% on average and was always higher than the negative trends, excepting for the ToC-BFAST (17%) map. Positive trends were more frequently found by the SWATI indicators (42% and 37% for AM and ESPI) and less frequently in LTT indicators (18% and 15%).

The difference between methods was always more important than the difference between using the NDVI annual means or ESPI time series, however, there are some small scale and intensity differences. The results also show that the Long-Term Trend method is more conservative, since it may fail to detect changes in vegetation productivity that occur in small portions of the time series or with a non-linear characteristic. The western montane steppes and grasslands with arid conditions showed the highest percentages of area with positive trends, whereas the Chaco dry forest and temperate grasslands and shrublands presented the highest percentages of area with negative trends, as shown in Table A2.

3.2. Expert Opinion Results

Remote sensing offers objective assessments of changes in biophysical properties measured from space. In most cases, datasets that can be used to assess the validity of those remotely sensed estimates are scarce, especially for validating changes over time. Using expert opinion, we were able to partly bridge that gap. A total of 131 expert opinions were submitted through the web application, covering all the phytogeographic provinces of Argentina. Of the six primary productivity maps presented, the LTT-AM was the most frequently selected as the most appropriate (27.5% of the total votes), followed by the SWATI-AM (22.9%) and the SSWATI-ESPI (16.8%). ToC-BFAST was the least selected map at a national scale.

At a regional level, as shown in

Table 2. Number of times each indicator was selected by experts in different ecoregions. Colors indicate the ranking of the indicator in each region: from blue (more votes) to red (less votes).

Biomes	LTT AM	LTT ESPI	SWATI AM	SWATI ESPI	ToC BAFST	SSWATI ESPI	Total
Temperate Grasslands, Savannas, and Shrublands (1,744,745 km2)	21	7	15	7	9	9	68
Montane Grasslands and Shrublands (272,555 km2)	3	1	2	0	1	6	13
Subtropical Moist Broadleaf Forests (112,711 km2)	7	5	2	2	0	5	21
Subtropical Grasslands, Savannas, and Shrublands (650,388 km2)	5	7	11	3	1	2	29
TOTAL (2,780,400 km2)	36	20	30	12	11	22	131

, LTT-AM was the most selected map both in temperate grasslands, savannas, and shrublands and in subtropical moist broadleaf forests. In subtropical grasslands, savannas, and shrublands, most experts chose the SWATI-AM and in montane grasslands and shrublands, the SWATI-ESPI was the most chosen.

Using the web-tool, experts also provided ancillary information about the main causes for the observed pattern. For the negative trends the most important were: deforestation, overgrazing, and overexploitation of the environment; while for the positive trends they were: recovery of burnt/disturbed areas, rainfall increase, and increase in protected areas. Trend analysis does not provide information on the causes, highlighting the importance of getting expert knowledge consultation on the indicators. Expert opinions found that some land changes associated with degradation processes may show positive trends, such as exotic afforestation of natural grasslands, invasive plant expansions, expansion of irrigation and alterations of water cycle, climatic-related events like snowmelt or shrinking water bodies and regrowth or greenness in fire-managed landscapes, native grassland to agriculture, and abandonment of degraded land to shrubland.

3.3. Detection of Negative Primary Productivity Trends in Plots with Forest Loss

The mean TI in areas which experienced tree cover loss varied significantly among methods. Also, all interaction terms among region, method, and type of land were significant in the adjusted model (p values < 0.001), showing that the alternative methods perform differently in every region and type of forest. However, in all regions and for all methods, areas with loss of forest lands showed significantly lower TI than in areas with loss of other woodlands, indicating that it is easier, as expected, to detect deforestation in areas with more canopy cover and higher trees. Since loss of forest lands and of other woodlands are a type of landcover and land-use change that should be reflected by a decrease in productivity, negative and lower values of TI were expected in these plots.

In general, average TI was negative and lower than −0.5, indicating light negative trends in deforested plots. However, in the subtropical moist forests (Paranaensis forests), average TI was higher than −0.5, and in some cases, it was even positive, indicating that deforestation in this region is hardly detected (mean TIs between −0.5 and 0.5 would correspond to a No Trend category). This could be due to the fact that deforestation in this region is performed to remove native tree species for planting fast growing exotic tree species for wood production. Differences among regions were also statistically significant, with subtropical Chaco dry forests showing the lowest levels of TI, as shown in Figure 4.

Average TI in all plots where deforestation was recorded was −0.82 for the SWATI-AM map, −0.78 for LTT-ESPI, and −0.78 for SWATI-ESPI (Figure 5). These three maps were regarded as the best performers in detecting negative trends of primary productivity in deforested plots, with SWATI-AM showing the lowest values. The LTT-AM average TI was −0.51, whereas SSWATI-ESPI and ToC-BFAST average TI were higher than −0.5. These results show that when calculating LTT indicators, the use of the ESPI time series instead of annual mean NDVI increases the ability to detect forest loss in all the regions.

When each type of forest and region were considered, there was no single best indicator to detect decreases in vegetation productivity due to deforestation, as shown in Figure 4. However, SWATI-AM presented significantly lower values of TI overall, reflecting a higher ability to detect loss of forest cover. When results were analyzed regionally, SWATI-AM was also the most sensitive indicator in subtropical dry forests (Chaco) and montane moist forests (Yungas). In temperate grasslands and shrublands (Espinal) and subtropical moist forests (Paranaensis forests), indicators using the ESPI time series performed better at detecting deforestation. For loss of forest lands, as shown in Figure 4a, the three best performers (SWATI-AM, SWATI-ESPI, and LTT-ESPI) showed similar behavior (lowest TIs). In the subtropical moist forests, where deforestation of native forests is frequently followed by the plantation of exotic tree species, the methods using time steps (SWATI and SSWATI) showed on average more ability to detect changes in deforested plots. In regions where land is cleared for agriculture or livestock rearing, all methods performed better. In the temperate shrublands of the Espinal region, there is a marked natural seasonality and the change to pasture or cultivated land is reflected more in the intra-annual variability than in the mean annual productivity and it was better captured when the ESPI time series were analyzed.

4. Discussion

To estimate the proportion of land that is degraded over total land area, it is necessary to take into account trends in land cover (SO1-1), trends in land productivity (SO1-2), and trends in carbon stocks (SO1-3), which increase the chances of detecting changes in land degradation. During the reporting process of the PRAIS 3 in Argentina, the trend in land productivity (SO1-2) was the most important indicator (accounting for 98% of the detected changes). It is expected that in large and diverse countries, the SO1-2 indicator would be the most reliable and easy to acquire indicator at a national scale, especially in comparison with the high cost and technical difficulty of obtaining multiple high quality land-cover and Soil Organic Carbon (SOC) maps. This is why, for the LDN baseline in Argentina, the changes in land productivity (SO1-2 indicator) were the focus of a validation process incorporating expert knowledge. However, estimating areas directly from a map may differ greatly from the true area because of map classification errors [42]. Therefore, it is recommended that area estimates are error-adjusted and accompanied by confidence intervals that quantify the sampling variability of the estimated area. To achieve this, a reference dataset with more reliability than the map must be obtained. There are several alternatives and good practice recommendations to design and implement area estimation methods with accuracy assessments to maps using remote sensing [43]. However, prior to these assessments, it is necessary to obtain a reliable base map in order to have meaningful results and low levels of uncertainties. Our results show that great differences arise among the various ways of computing trends in primary productivity, highlighting the importance of integrating expert opinions to choose the more adequate maps at national and local levels before obtaining reference sample data for the estimation of areas and their uncertainties.

The LTT-AM was the map most frequently (27%) identified as the one that better reflected changes in land productivity. This method is the most widely used trajectory approach and is implemented in Trends.Earth [20], a free and open source tool that uses earth observation systems to evaluate changes in land condition and that can be used to a national level reporting to SDG Indicator 15.3.1. The LTT-AM was followed by the SWATI-AM map, which was chosen 22.9% of times. The rest of the methods were selected less than 16.8% of times. The SWATI approach, presented in this paper, was designed to capture both short-term and long-term trends, in contrast to LTT-AM, which may fail to detect changes with short-term persistence. Analyzing data focusing only on short-term or long-term persistence may affect the trend-detection performance [44]. In long periods of time, linear trends tend to get smaller and non-linear behavior trends will also tend to be not significant. That is probably the reason why the LTT indicators show less diversity of behaviors than the methods which consider the short-term persistence. The SWATI-AM map showed more variability and had a good performance on both hard and soft validation datasets, particularly, by detecting decreases of primary productivity due to forest loss.

It is expected that not all the land-cover/land-use changes will be detected by each trend detection method, as shown by global analysis on the LTT-AM and World Overview of Conservation Approaches and Technologies (WOCAT) SLM database in Gonzalez-Roglich et al. [22]. Our results from the analysis on the deforested plots confirm this, as the best indicator varied among ecoregions. Results from expert perception may always include some degree of bias at this scale [6], but when added in, the agreement and disagreement among the alternative maps help to better understand land processes and to narrow down the more adequate maps. A multi-indicator approach including expert knowledge interpretation is probably the best way to approximate a more reliable map of productivity trends. Future studies could advance in combining different methods, perhaps with an ensemble modeling approach [45,46,47], to try to capture the best features of each indicator into a single product. Also, comparisons and inclusion of other different approaches easy to implement in cloud computing environments, like BACI (Before-After, Control-Impact) [48] or Trends.Earth [20], could prove very beneficial for these types of objectives.

At a regional scale, in some ecosystems, and especially on those cases where land changes have impact on seasonality but do not impact the overall productivity, the use of the ESPI time series was favored over the use of AM, as shown in Appendix A. This type of behavior was also suggested in the work of de Jong et al. [49] and Paruelo et al. [28]. When the LTT model was considered, the use of ESPI over the AM time series created more diverse responses, as shown in Figure 2. In these cases, to just analyze the change in net annual productivity may not be enough and intra-annual variations should be considered to better assess the potential degradation or improvements. Regarding the least selected map, it is important to consider that the ToC-BFAST method was originally obtained using the 16-day dataset, including higher intra-annual information, which produced different metrics and results, but was later forced into the trend intensity (TI) categorization to make it comparable with the other maps in the expert consultation. In that conversion some important information was lost, like the year of the break point in the trend, which may provide important information for local and sub-national level analysis. Also, the type of change class was not used in this analysis, which is an additional option to present high frequency data at a national scale rather than the trend intensity.

On average, the six alternative methods detected decreases of land productivity in 20% of the Argentine territory. However, according to the LTT-AM, which was the most selected method, 9% of the territory presented negative primary productivity trends, most of which (6.73%) correspond to decreases of less than 25% of the initial value (light negative trends). According to the SWATI-AM map, which was also highly selected by the experts, 23.5% of the national territory presented a negative trend, from which 14% correspond to light negative trends and 9% to moderate negative trends. The high variability among percentages of land with negative trends found among methods (from 9% to 35%) has a direct impact on the construction of Indicator 15.3.1 in the context of UNCCD reporting for the 2030 Agenda, as well as for Argentina’s LDN monitoring. The difference regarding percentage of land with negative trends between the two most selected maps (LTT-AM and SWATI-AM) is 14.5% of Argentine territory, which corresponds to approximately 400,000 km², more than Germany´s area. Nevertheless, the validation methods aided in decreasing this range, since they allowed us to obtain a more solid set of indicators and a lower level of associated uncertainty.

In Argentina, previous reports indicated 38% of land with negative trends [26,50] for the 2000–2014 time interval using a very similar approach to LTT-AM. Our results show that using the same LTT-AM approach on the 2000–2018 period, the area with negative trends dropped to 9%. In monitoring the country´s progress towards LDN, this could be translated as an improvement in the country combating desertification or it may be related to the LTT method´s sensitivity to the choosing of the start/end period and longitude of the series [51]. Future studies of the method’s sensitivity to different time intervals could be very beneficial in finding good indicators for LDN monitoring.

By adding the trend intensity categorization based on trend slopes (magnitude of change by unit of time), we were able to provide information on how the different indicators not only detect, but also characterize, the time series behavior. Normally, the SO1-2 indicator is reported in terms of three categories (negative, positive, or neutral), since it is part of SDG Indicator 15.3.1, where the one out all out method [9] is applied, and negative trends are directly associated with degraded land. From a country or sub-national perspective, to identify if the trend is strong or light may provide stakeholders a chance to better allocate resources and priorities. In the monitoring phase, this information may also help to identify acceleration or deceleration of the trend intensity, i.e., a change from strong negative to light negative. This may be important for decision makers, especially when SLM practices or interventions are being evaluated. In this work, the slope was categorized as relative to the starting value, in order to have a comparable measure across different ecosystems and climatic zones. This method could potentially be improved to incorporate some elements of the net change in absolute values. For example, increasing the NDVI value from 0.08 to 0.12 is 50% and a 0.04 NDVI unit increase, while going from 0.6 to 0.9 is also 50% but is actually 0.3 NDVI units. Both changes may be important for their own environment by using relative change (percentage) values, but when the absolute change is considered, the amount of productivity and the associated carbon that is potentially fixed, is absolutely different in both scenarios. Future studies would be recommended to establish a proper way to make categories that could capture both the relative and absolute change to provide further information for the decision making process.

Computational effort is also a main issue that sometimes becomes a limiting factor when calculating indicators over large areas. To solve this issue, the LTT, SWATI, and SSWATI indicators were developed for this study to be computed for any country or region of the world using the Google Earth Engine platform, which has a demonstrated power of analysis [52]. In the Supplementary Material, the provided script will allow other users to obtain the primary productivity trend maps in a fast and inexpensive way for their study area.

5. Conclusions

As different maps captured different land processes, including land-use/land-cover changes, expert and ground knowledge needs to be considered in order to choose the best indicator of primary productivity trend and to evaluate whether the trend indicates actual land degradation or land improvement. The tool developed to make the expert consultation, built on open source software, was reliable, transparent, and secure, and was simple and cheap to implement at a national scale. It allowed us to collect expert knowledge in a systematic and objective way, facilitating the analysis of the collected information. A substantial part of the success in combating land degradation and achieving LDN is about finding consensus between different stakeholders. The provision of alternative methods and tools that enable sharing opinions and participatory discussions, improved the communication and baseline assessment during the LDN target setting in Argentina. Further studies that promote discussion on different methods to obtain better indicators could increase the chances of achieving LDN and improve planning for 2030 Agenda activities. Our tools and processes resulted in more people participating, better discussions, deeper understanding of the challenges ahead and at the end, more appropriation of the whole participatory process, which is key to building consensus.