A Multiscale Productivity Assessment of High Andean Peatlands across the Chilean Altiplano Using 31 Years of Landsat Imagery

: The high Andean peatlands, locally known as “bofedales”, are a unique type of wetland distributed across the high-elevation South American Altiplano plateau. This extensive peatland network stores signiﬁcant amounts of carbon, regulates local and regional hydrological cycles, supports habitats for a variety of plant and animal species, and has provided critical water and forage resources for the livestock of the indigenous Aymara communities for thousands of years. Nevertheless, little is known about the productivity dynamics of the high Andean peatlands, particularly in the drier western Altiplano region bordering the Atacama desert. Here, we provide the ﬁrst digital peatland inventory and multiscale productivity assessment for the entire western Altiplano (63,705 km 2 ) using 31 years of Landsat data (about 9000 scenes) and a non-parametric approach for estimating phenological metrics. We identiﬁed 5665 peatland units, covering an area of 510 km 2 , and evaluated the spatiotemporal productivity patterns at the regional, peatland polygon, and individual pixel scales. The regional assessment shows that the peatland areas and peatlands with higher productivity are concentrated towards the northern part of our study region, which is consistent with the Altiplano north–south aridity gradient. Regional patterns further reveal that the last seven years (2011–2017) have been the most productive period over the past three decades. While individual pixels show contrasting patterns of reductions and gains in local productivity during the most recent time period, most of the study area has experienced increases in annual productivity, supporting the regional results. Our novel database can be used not only to explore future research questions related to the social, biological, and hydrological inﬂuences on peatland productivity patterns, but also to provide technical support for the sustainable development of livestock practices and conservation and water management policy in the Altiplano region.


Introduction
Although wetlands cover only about 6% of the global land surface [1,2], they support remarkable levels of biological diversity and provide a wide spectrum of ecological goods and services [3]. Among the diverse types of wetlands, peatlands stand out as extraordinary carbon sinks with estimated carbon stocks between 113 and 612 Pg (billion tonnes) globally [4,5]. While peatlands are primarily   Among the diverse types of wetlands, peatlands stand out as extraordinary carbon sinks with estimated carbon stocks between 113 and 612 Pg (billion tonnes) globally [4,5]. While peatlands are primarily found in the Northern Hemisphere boreal and subarctic regions [5], a unique peatland system, sustained by groundwater flows and snow-and glacier-melt, spans across the Central Andes Altiplano (~14°S-26°S)-the largest high-elevation, semiarid plateau in South America (>4000 m.a.s.l.) (Figures 1 and 2). Due to limited water availability, indigenous communities have depended on the extensive high Andean peatland network, locally known as "bofedales", as water reservoirs and livestock forage resources for thousands of years [6][7][8].   The high Andean peatlands have the highest primary productivity among the distinct terrestrial ecosystems of the Altiplano region, with average carbon accumulation rates at least an order of magnitude higher than Northern Hemisphere peatlands [9][10][11]. Many of these peatlands represent valuable paleoclimate archives, as they have accumulated peat for up to 10,000 years and can reach depths of up to 10 meters [11][12][13][14]. Since climate and land use change has resulted in the degradation of these fragile ecosystems, and because peatlands have the potential to shift from carbon sinks to sources, there is a great need to evaluate changes in peatland extent and productivity [4]. This is especially true in semiarid regions such as the Central Andes Altiplano, where peatlands play important hydrological and ecosystemic roles [15,16] and are highly vulnerable to the projected changes in the regional hydroclimate [17]. The demonstrated capability to map these peatlands in detail [18], and to use satellite-based vegetation indices as indicators of their primary productivity [19][20][21], allows us to develop regional-scale studies of the Altiplano peatlands utilizing a fine-scale spatiotemporal remote sensing approach.
A number of studies at different spatial and temporal scales have used satellite-derived normalized difference vegetation index (NDVI) data, along with other vegetation indices, to study different aspects of the high Andean peatlands. These studies include, but are not limited to, evaluations of spatiotemporal changes in primary productivity [19][20][21][22], classifications of migratory and endemic bird habitats [23], and analyses of spatiotemporal eco-hydrological patterns [19,24,25]. For instance, Moreau et al. [19] assessed the temporal biomass dynamics of the Andean peatlands using NOAA/AVHRR NDVI data and found a close relationship between NDVI and humid (green) total peatland biomass measured at four Northern Bolivian Altiplano sampling sites. Otto et al. [24] developed a hydrological peatland classification scheme using 30 m Landsat NDVI and the normalized difference infrared index (NDII) to distinguish between perennial and temporal subtypes of peatlands in a bounded area of the Northern Altiplano in Peru. Using 250 m MODIS Terra (MOD13Q1) data, Otto et al. [24] also explored temporal changes in productivity. Izquierdo et al. [24] classified Landsat TM scenes from 2009 to 2013 in the Southeastern Altiplano in Argentina to map peatland density, size, and altitude, and to classify sub-watersheds throughout the region. Baldassini et al. [22] focused more broadly on characterizing distinct vegetation types, including peatlands, throughout the Southeastern Altiplano in Argentina to assess differences in aboveground net primary productivity. Similar to Otto et al. [24], Baldassini et al. [22] used data from MODIS (MOD13Q1) and was limited to a ten-year analysis period (2000-2009). Dangles et al. [21] built on these studies and analyzed a combination of Landsat imagery from 1984 to2011 and PLEIADES images from 2013 to monitor the number and size of peatlands in the glaciated Cordillera Real in the Northeastern Bolivian Altiplano. Likely caused by the increased melting of glaciers, Dangles et al. [21] detected an overall increase in peatland cover over time. While these studies demonstrate the ability to leverage remote sensing time series to assess and monitor the high Andean peatlands, they also highlight the importance of encompassing greater spatial and finer temporal resolutions to evaluate the ecosystem responses to environmental changes across the Central Andes. At present, studies surrounding the high Andean peatlands located within the vast western region of the Altiplano are scarce. This region is located in the Chilean Andes, adjacent to the Atacama Desert, and spans more than 60,000 km 2 , representing the driest area of the Altiplano plateau.
A better understanding of peatland dynamics in this region is critical due to current and projected environmental changes. Instrumental registries already show marked changes in the Altiplano climate, including less frequent rainfall events [26] and a positive warming trend over the last 50 years [27]. Furthermore, climate models project future increases in aridity for the region, especially across the western Altiplano in Chile where changes would be most severe [17,28]. In this regard, it is of particular interest to study the productivity dynamics and distribution of the dry western Altiplano peatlands. In the present study, we develop the first large-scale and high-resolution multidecadal database for peatland monitoring across the western Altiplano in Chile (63,705 km 2 ). By compiling 31 years of Landsat NDVI records, we: 1) develop a detailed inventory of active (green) high Andean peatlands in this semiarid region, 2) characterize the geographical patterns in productivity at different spatial scales Remote Sens. 2019, 11, 2955 4 of 23 across an~800 km latitudinal gradient along the Andes, and 3) assess the regional spatiotemporal changes in productivity using time series analysis and phenological reconstructions. We expect that the database and findings presented here will substantially improve our understanding of the spatiotemporal patterns in Andean peatland productivity during the last three decades and can be used for more detailed peatland assessments across the study region in the future.

Satellite Data
We used all available Landsat Surface Reflectance data (Collection 1 Level 2) from the 5-TM, 7-ETM+, and 8-OLI to prepare an NDVI raster stack for the entire study area. A total of 8997 scenes were requested and downloaded from the USGS Earth Explorer portal, corresponding to 10 Path/Rows (002/72:73, 001/72:75, 233/74:75 and 232/76). The quality assessment bands of these products were used to discard pixels flagged as snow/ice, cloud, cloud shadow, and cirrus detected with high and medium confidence (see Landsat 8 and Landsat 4-7 Surface Reflectance Product Guide, available at the USGS Earth Explorer portal https://earthexplorer.usgs.gov/). Of the 56.3 billion pixels from our Landsat subset, 10.5 billion (19%) were discarded using these quality assessment bands. The resulting NDVI time series raster stack was comprised of 31 Southern Hemisphere growing seasons (GSs) from 1986-87 to 2016-17.

Study Region
The Altiplano plateau extends over 1000 km latitudinally (15 • S-23 • S) and is situated in the Central Andes between the Atacama Desert to the west and the tropical humid lowlands to the east. In this semiarid bioregion, annual precipitation decreases from around 450 millimeters (mm) in the northeast region bordering the Amazon basin to less than 150 mm in the southwest sector adjacent to the Atacama Desert. More than 80% of the total annual precipitation falls during summer as part of the South American monsoon [29], which is partially modulated by the El Niño-Southern Oscillation phenomenon [30]. Our study region encompassed the entire Chilean Altiplano, representing the driest sector of the plateau, and a 30 km buffer beyond the country border to prevent the arbitrary division of peatlands that cross administrative country boundaries ( Figure 1). In total, our study covered an area of 63,705 km 2 , above 3700 m.a.s.l., and between 17.3 • S and 24.0 • S and 70.0 • W and 66.8 • W. Floristically, the high Andean peatlands are primarily composed of cushion-forming species from the Juncaceae family such as Distichia muscoides, Oxychloe andina, and Patosia clandestine [31]. Measurements within our study area made by Cooper et al. [10] found that peatland height increases from 0.96 to 5.37 cm/yr, has a mean bulk density of 0.081 g/cm 3 , and mean organic carbon production rates ranging from 1.5 to 4.0 kg/m 2 /yr.

Preliminary Extent Assessment of Six Pilot Peatlands
The temporal productivity dynamics of the high Andean peatlands have been evaluated using either shifts in area and peatland extent [18,21] or changes in green biomass within a fixed area [19]. For the second approach, defining the appropriate fixed peatland area is critical. The maximum historical extent is often preferred as it captures the NDVI differences related to contractions of the peatland area and to reductions of green biomass within the specified area. Here, we use the fixed area approach based on the maximum historical extent. Prior to executing our classification procedure for the entire study region, we carried out a preliminary assessment of the monthly variation in extent for six pilot peatlands ( Figure 1). All of the pilot peatlands showed seasonal changes in area (Figure 3), reaching a maximum extent between February and March (towards the end of the Southern Hemisphere summer). Years with the highest peatland green extent were 1989,1990,1991,1996,2000,2001,2002,2004,2005,2006 and 2015 (see Figure 3, where the area time series of three peatlands are depicted). For this preliminary exercise, we analyzed cloud-free Landsat data from 1986 to 2017 for the following 6 pilot peatlands: from north to south, Cosapilla, Pomerape, Guallatire, Lirima, Tatio, and Machuca ( Figure 1). A manual delineation of the pilot peatlands was carried out using Google Earth for two seasons (Feb/March and Sep/Oct) in both 2015 and 2016. This was cross-checked with the field observations and the photographs obtained during the field campaigns. We then extracted the NDVI values of all Landsat pixels inside each of the manually delineated peatland polygons from the nearest available dates corresponding to the 2015 and 2016 spring and winter seasons. For each peatland, we used boxplots to assess the distribution of the NDVI values. Based on these boxplots, we set an NDVI threshold to automate the delineation of the pilot peatlands for each available Landsat scene. Using these time series, we determined the maximum historical extent of each pilot peatland between 1986 and 2017. For each peatland, we used boxplots to assess the distribution of the NDVI values. Based on these boxplots, we set an NDVI threshold to automate the delineation of the pilot peatlands for each available Landsat scene. Using these time series, we determined the maximum historical extent of each pilot peatland between 1986 and 2017.

Large-scale Digital Inventory of High Andean Peatlands
Building on the results from the pilot peatland assessment, we delineated all peatlands across the study region. To determine the maximum extent for each individual peatland (segmentation process), we used 50-55 Landsat scenes, depending on the availability for each path/row, from the years in which the pilot peatlands showed the greatest spatial extent ( Figure 3). Figure 4 provides an example of the full segmentation procedure for the Guallatire peatland. Following a Boolean reclassification of each of the 50-55 scenes (1 for "peatland" and 0 for "non-peatland") based on a simple thresholding procedure, we delineated the peatland polygons for the digital inventory using pixels classified as peatlands in at least 75% of the 50-55 scenes (Figure 3c-d). The 75% threshold was selected to account for missing data due to cloud cover. Polygons comprised of fewer than five pixels (4500 m 2 ) were excluded from our digital inventory. Finally, an object-based accuracy assessment was done by comparing manually delineated peatlands (using Google Earth imagery) and automatically segmented peatlands (using Landsat imagery, this study) [32]. For 1% of all segmented peatlands, a total delineation error (ranging from 0 to 1) was calculated, which accounts for both overestimated and underestimated area [32]. We compared the Landsat-segmented peatlands with the manually delineated peatlands at dates when the available high-resolution images from Google Earth were closest in time to the specific Landsat scenes. Since the spatial resolution of Google Earth imagery allows for a more detailed delineation of peatlands than the Landsat spatial resolution (30

Large-Scale Digital Inventory of High Andean Peatlands
Building on the results from the pilot peatland assessment, we delineated all peatlands across the study region. To determine the maximum extent for each individual peatland (segmentation process), we used 50-55 Landsat scenes, depending on the availability for each path/row, from the years in which the pilot peatlands showed the greatest spatial extent ( Figure 3). Figure 4 provides an example of the full segmentation procedure for the Guallatire peatland. Following a Boolean reclassification of each of the 50-55 scenes (1 for "peatland" and 0 for "non-peatland") based on a simple thresholding procedure, we delineated the peatland polygons for the digital inventory using pixels classified as peatlands in at least 75% of the 50-55 scenes (Figure 3c,d). The 75% threshold was selected to account for missing data due to cloud cover. Polygons comprised of fewer than five pixels (4500 m 2 ) were excluded from our digital inventory. Finally, an object-based accuracy assessment was done by comparing manually delineated peatlands (using Google Earth imagery) and automatically segmented peatlands (using Landsat imagery, this study) [32]. For 1% of all segmented peatlands, a total delineation error (ranging from 0 to 1) was calculated, which accounts for both overestimated and underestimated area [32]. We compared the Landsat-segmented peatlands with the manually delineated peatlands at dates when the available high-resolution images from Google Earth were closest in time to the specific Landsat scenes. Since the spatial resolution of Google Earth imagery allows for a more detailed delineation of peatlands than the Landsat spatial resolution (30 m), we resampled the Google Earth-delineated polygons to Landsat resolution prior to calculating the delineation errors.

Time-series Analysis and Productivity Assessment at Different Spatial Scales
We evaluated NDVI time series at the regional, individual peatland, and pixel scales using our full digital inventory. The regional level represents the average of all peatland pixels across the study region. Individual peatland time series are calculated as the average of all pixels within each polygon. Pixel analyses include the NDVI values for all 30 × 30 meter pixels classified as peatlands. For each spatial scale, we reconstructed the annual phenological cycle using a non-parametric approach implemented in the "npphen" R package [33]. This tool has proven to be effective for studying vegetation phenological metrics using large remote sensing datasets [34][35][36][37]. With npphen, all NDVI observations (Figure 5a) are arranged by the phenological year (from July to the following June for the Southern Hemisphere) (Figure 5b). This so-called NDVI-day of growing season (DGS) space is used to generate a kernel density estimation for NDVI at all DGSs (Figure 5c). We estimated the annual vegetation productivity by calculating the area under the most likely observed NDVI value curve, a phenometric commonly used as a productivity proxy [19,38,39]. We use this accumulated annual NDVI sum (NDVISUM) to study peatland productivity across the western Altiplano at the pixel, polygon, and regional scales.

Time-Series Analysis and Productivity Assessment at Different Spatial Scales
We evaluated NDVI time series at the regional, individual peatland, and pixel scales using our full digital inventory. The regional level represents the average of all peatland pixels across the study region. Individual peatland time series are calculated as the average of all pixels within each polygon. Pixel analyses include the NDVI values for all 30 × 30 meter pixels classified as peatlands. For each spatial scale, we reconstructed the annual phenological cycle using a non-parametric approach implemented in the "npphen" R package [33]. This tool has proven to be effective for studying vegetation phenological metrics using large remote sensing datasets [34][35][36][37]. With npphen, all NDVI observations (Figure 5a) are arranged by the phenological year (from July to the following June for the Southern Hemisphere) (Figure 5b). This so-called NDVI-day of growing season (DGS) space is used to generate a kernel density estimation for NDVI at all DGSs (Figure 5c). We estimated the annual vegetation productivity by calculating the area under the most likely observed NDVI value curve, a phenometric commonly used as a productivity proxy [19,38,39]. We use this accumulated annual NDVI sum (NDVI SUM ) to study peatland productivity across the western Altiplano at the pixel, polygon, and regional scales.

Spatiotemporal Changes in Productivity
We first calculated the regional monthly mean NDVI time series based on all peatland pixels in our study area. Then, we removed the seasonal component from the regional NDVI signal using seasonal-trend decomposition implemented in the DBEST R package, which is based on the LOESS (locally estimated scatterplot smoothing) method [40], a filtering procedure that employs locally weighted regression for smoothing [41]. Finally, we assessed the significant regime shifts in the regional NDVI signal using the Rodionov method, which identifies abrupt changes from one relatively stable state to another through a sequential algorithm of t-tests [42]. For each regional regime shift, we analyzed the NDVISUM changes at the polygon and pixel levels to explore local geographical patterns within and between periods. We summarized the results of this highresolution analysis for the entire study region using productivity histograms based on 0.1° latitude bands.

Digital Inventory of High Andean Peatlands
Based on the boxplots from the manual delineation ( Figure 6), we used a threshold of NDVI >0.23 to automatically delineate the peatland polygons. Since the Altiplano is a cold, semi-arid region, only peatland vegetation registers NDVI values greater than 0.23, justifying the simple thresholding procedure (Figures 4 and 6).

Spatiotemporal Changes in Productivity
We first calculated the regional monthly mean NDVI time series based on all peatland pixels in our study area. Then, we removed the seasonal component from the regional NDVI signal using seasonal-trend decomposition implemented in the DBEST R package, which is based on the LOESS (locally estimated scatterplot smoothing) method [40], a filtering procedure that employs locally weighted regression for smoothing [41]. Finally, we assessed the significant regime shifts in the regional NDVI signal using the Rodionov method, which identifies abrupt changes from one relatively stable state to another through a sequential algorithm of t-tests [42]. For each regional regime shift, we analyzed the NDVI SUM changes at the polygon and pixel levels to explore local geographical patterns within and between periods. We summarized the results of this high-resolution analysis for the entire study region using productivity histograms based on 0.1 • latitude bands.

Digital Inventory of High Andean Peatlands
Based on the boxplots from the manual delineation ( Figure 6), we used a threshold of NDVI >0.23 to automatically delineate the peatland polygons. Since the Altiplano is a cold, semi-arid region, only peatland vegetation registers NDVI values greater than 0.23, justifying the simple thresholding procedure (Figures 4 and 6). We identified a total of 5665 peatland polygons distributed between 17.3°S and 23.7°S, covering a total area of 510.27 km 2 (Figures 7 and 8). Detailed maps of all identified peatlands in our study area are provided in the Appendix. Of the total peatland area, 73% is found in the more humid northern latitudes between 17.3°S and 19.0°S, where the annual rainfall can reach up to ~450 mm. The peatland area significantly decreases south of 19.0°S and nearly disappears around 21°S. A second and smaller cluster of peatlands is located toward the much more arid southern margin of our study region (22.5°S-23.7°S). We mapped 3325 peatland polygons within the Chilean administrative borders and an additional 26 units that cross into the surrounding countries, summing to a total area of 291.6 km 2 . The mean size of all peatland polygons was 0.09 km 2 (about 100 Landsat pixels) with a standard deviation of 0.0045 km 2 . The largest peatland in our inventory is found at the bases of the Bolivian Sajama and Parinacota volcanoes with an area of 33.66 km 2 . This is followed by Pomerape in Chile (19.83 km 2 ) and Nasahuento on the Chilean-Bolivian border (15.41 km 2 ) ( Figure 9). In general, the largest peatlands form at the base of extensive valleys. Roughly 70% of the peatland area is located between 4000 and 4500 m in elevation and almost 80% of the peatland area is found in areas with null or low slope (0-5%), following the direction and shape of the drainage network. Larger and more productive peatlands are concentrated in the northern cluster, while smaller, less abundant, and less productive peatlands are found in the southern cluster (Figures 8 and 10). We identified a total of 5665 peatland polygons distributed between 17.3 • S and 23.7 • S, covering a total area of 510.27 km 2 (Figures 7 and 8). Detailed maps of all identified peatlands in our study area are provided in the Appendix A. Of the total peatland area, 73% is found in the more humid northern latitudes between 17.3 • S and 19.0 • S, where the annual rainfall can reach up to~450 mm. The peatland area significantly decreases south of 19.0 • S and nearly disappears around 21 • S. A second and smaller cluster of peatlands is located toward the much more arid southern margin of our study region (22.5 • S-23.7 • S). We mapped 3325 peatland polygons within the Chilean administrative borders and an additional 26 units that cross into the surrounding countries, summing to a total area of 291.6 km 2 . The mean size of all peatland polygons was 0.09 km 2 (about 100 Landsat pixels) with a standard deviation of 0.0045 km 2 . The largest peatland in our inventory is found at the bases of the Bolivian Sajama and Parinacota volcanoes with an area of 33.66 km 2 . This is followed by Pomerape in Chile (19.83 km 2 ) and Nasahuento on the Chilean-Bolivian border (15.41 km 2 ) (Figure 9). In general, the largest peatlands form at the base of extensive valleys. Roughly 70% of the peatland area is located between 4000 and 4500 m in elevation and almost 80% of the peatland area is found in areas with null or low slope (0-5%), following the direction and shape of the drainage network. Larger and more productive peatlands are concentrated in the northern cluster, while smaller, less abundant, and less productive peatlands are found in the southern cluster (Figures 8 and 10). Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 24

Accuracy Assessment of the Digital Inventory
The object-based accuracy assessment showed a mean delineation agreement of 0.79 (standard deviation of 0.1) between the Landsat-delineated peatlands and the manually delineated peatlands using Google Earth imagery. This metric, ranging from 0 to 1, accounts for both the overestimated and underestimated area for each peatland polygon in relation to the reference [32]. Larger disagreements were more common among the smaller peatlands, where fewer missing pixels cause larger differences in extent. In some cases, the inability to differentiate active peatlands from depleted and inactive peatlands on Google Earth imagery, without the clear contrast from the NDVI (see Figure 4), could also explain some of these disagreements.

Multiscale Annual Peatland Productivity Assessment
The reconstructed regional NDVI annual phenological cycle of the western Altiplano ( Figure 10) shows recurrent seasonal behavior and a mean productivity of 3.96 in terms of the NDVI SUM . Based on the kernel density estimation, the transition towards the growing season onset is more gradual and less variable than the growing season offset (Figure 10). The growing season begins around the day of the growing season (DGS) 94 (early October), reaching a maximum value around DGS 237 (late February), and maintaining high productivity through to DGS 258 (mid-March), which corresponds to the late summer season. Following DGS 258, NDVI drops rapidly until DGS 318 (mid-May), when the growing season ends. Standard to phenological studies, we define the growing season onset and offset as the dates when 20% of peak NDVI is reached [39,43]. At the pixel scale, annual peatland productivity reveals significant spatial heterogeneity within individual peatlands, which is likely a response to local hydrological conditions and/or livestock management. Higher NDVI SUM values are found towards the northern limits of our study area (Figure 9a), primarily influenced by the concentration of larger peatlands (Figure 9b,c). For example, the Pomerape and Sajama peatlands are responsible for the peak in productivity observed around 18.1 • S (Figure 9a).
The regime shifts analysis of the regional mean NDVI time series (Figure 11) indicates that Andean peatlands have had relatively stable productivity between 1986 and2011 with a mean NDVI value of roughly 0.34 (Figure 11c). However, this is interrupted by two drops and recoveries in NDVI of a magnitude of about 0.04 in 1989-1993 and 1998-2000. In June 2011, we detected another regime shift with a positive NDVI magnitude of 0.03. This green-up has been maintained through to 2017 and represents the period with the highest productivity between 1986 and 2017. For each of the six distinct periods identified by the regime shift analysis of the regional NDVI time series, we calculated the NDVI SUM at the pixel level and summarized the relative peatland area (%) in five productivity (NDVI SUM ) classes ( Figure 12). In the sixth period (2011-2017), 17% of the total peatland area falls within the highest productivity class (NDVI SUM > 6), reaching a maximum extent in terms of area for this productivity class. This is followed by the first and second periods (13% each). For periods with the least productivity (2 and 4), less than 9% of the total area is included in the highest productivity class. While the percentage of area with the highest annual productivity (NDVI SUM > 6) in the sixth period increased by 25% with respect to the first period, the area with the lowest annual productivity (NDVI SUM < 3) decreased by 18% ( Figure 11). Overall, the extreme classes of productivity have changed the most drastically since the first period. However, more areas have witnessed gains in annual productivity than losses, which explains the recent regional greening detected by the regime shift test. middle section show negligible changes for most periods. This multiscale analysis demonstrates how our digital inventory can be used not only to study regional changes in productivity, but also patterns within individual peatlands. While we cannot refer to all peatlands individually in this study, we believe that our inventory is a powerful tool for the scientific community, stakeholders, and local managers for further productivity analyses.   When disaggregating these productivity differences latitudinally (Figure 13), it is clear that the two periods with drops in productivity (P2 and P4) are dominated by areas with negative productivity differences across all latitudes with respect to P1. P4 shows the largest area of negative productivity differences, despite 10% to 25% of the area, considering different latitudinal ranges, presenting positive changes. In the most recent period (P6), peatlands show contrasting productivity differences, with both negative and positive differences for all latitudinal bands. Aside from a narrow band between 22.2 • S and 22.3 • S, there is a higher percentage of area with positive differences, evidencing overall greening. The presence of both negative and positive productivity differences indicates that despite this general pattern, there is heterogeneous variability at the local scale. In order to illustrate this local spatiotemporal heterogeneity, Figure 14 illustrates productivity differences at the pixel scale between all six periods for the Guallatire peatland. In general, we observe positive productivity differences in the southern part of the Guallatire peatland, while negative differences are more frequent in the central and northern areas near the basin headwaters. Some branches in the middle section show negligible changes for most periods. This multiscale analysis demonstrates how our digital inventory can be used not only to study regional changes in productivity, but also patterns within individual peatlands. While we cannot refer to all peatlands individually in this study, we believe that our inventory is a powerful tool for the scientific community, stakeholders, and local managers for further productivity analyses.

Discussion
Here, we provide the first digital peatland inventory for the entire western Altiplano (63,705 km 2 ; 17.3°S-24°S) and evaluate productivity at the regional, peatland polygon, and individual pixel scales. Using 31 years of Landsat data, this study is the most detailed spatiotemporal analysis of this ecosystem type within the Altiplano region. Since the Atacama Desert and the western Chilean Altiplano are among the regions with the lowest cloud fractions globally [44], conditions are advantageous for studying spatiotemporal patterns of peatland vegetation with the 30+ year catalogue of Landsat images. The medium spatial resolution of Landsat (30 m) and the significant amount of cloud-free data make this dataset particularly useful to study the high Andean peatlands,

Discussion
Here, we provide the first digital peatland inventory for the entire western Altiplano (63,705 km 2 ; 17.3 • S-24 • S) and evaluate productivity at the regional, peatland polygon, and individual pixel scales. Using 31 years of Landsat data, this study is the most detailed spatiotemporal analysis of this ecosystem type within the Altiplano region. Since the Atacama Desert and the western Chilean Altiplano are among the regions with the lowest cloud fractions globally [44], conditions are advantageous for studying spatiotemporal patterns of peatland vegetation with the 30+ year catalogue of Landsat images. The medium spatial resolution of Landsat (30 m) and the significant amount of cloud-free data make this dataset particularly useful to study the high Andean peatlands, especially considering that most peatlands are small (9 ha on average) and elongated (Figure 9 and Appendix A). For this reason, most remote sensing studies on high Andean peatlands have used Landsat data (i.e., [18,20,21,24]). Coarser data, such as the 250 m MODIS vegetation indices product, has been widely used for peatland mapping in the Northern Hemisphere where peatlands are larger in size [4]. However, this MODIS product has shown limitations and systematic biases due to spatial constraints for studying high Andean peatlands, even when downscaled [24]. Casagranda et al. [19] also found that MODIS constrained analyses of high Andean peatlands and had to use an algorithm to estimate the spatial contribution of peatlands to the NDVI values. In the present study, the Landsat catalogue allowed us to successfully identify the long-term maximum extent of peatland units and to assess productivity at distinct spatial scales. Because no other dataset can provide such a long record at this detailed spatial scale, the Landsat catalogue is a unique and valuable dataset for understanding the spatiotemporal dynamics of western Altiplano peatlands over the past 30 years.
In addition to updating our high Andean peatland inventory with new Landsat data, current monitoring satellites such as the Sentinel 2a and 2b of the Copernicus Programme of the European Space Agency can be integrated to provide a productivity assessment at an even finer spatial scale in recent years (2015-present) [45]. For the management of specific peatland units, the Sentinel 2 10-meter spatial and 5-day temporal resolutions support more detailed analyses of peatland productivity compared to Landsat, allowing local organizations to better evaluate grazing schemes, water management, and conservation practices at a much smaller scale. Commercial very high spatial resolution sensors (<5 m), including WorldView2, WorldView3, and GeoEye, among others, are not as appealing because they do not provide data on a regular basis. The RapidEye constellation may be an exception, with daily data acquisition and 5 m spatial resolution, but at a high cost. As a result, we believe that the combination of 30 m Landsat data and 10 m Sentinel 2 data are the most promising remote sensing sources for studying long-term high Andean peatland productivity (30+ years) at high spatiotemporal resolutions moving forward. These decadal studies can be further complemented with detailed assessments of peatland biomass [46,47], functional traits and plant composition [48,49], soil moisture [50,51] and water table measurements [52][53][54][55] using active radar or LiDAR systems, hyperspectral sensors, optical sensors on unmanned vehicles or multi-sensor approaches recently developed for other regions around the world.
Our productivity results show that peatlands with higher productivity are concentrated toward the northern part of our study area. This is consistent with the Altiplano north-south aridity gradient, where roughly double the amount of annual precipitation falls in the northern sector in comparison to the southern extreme. In addition to the regional gradient, there is high spatial heterogeneity in productivity within peatland units ( Figure 14). These variations could be related to differences in the sensitivity of vegetation at specific micro-sites stemming from unique interactions between drainage conditions, microtopography, and exposure. Individual peatlands are also modulated by the grazing strategies of indigenous communities, which could cause highly local spatiotemporal productivity variations within a peatland. However, a marked cyclical winter-summer pattern in productivity is observed for high Andean peatlands across all spatial scales, which is likely a response to the large-scale precipitation and/or temperature controls on the Altiplano region. Considering the complex dynamics of these ecosystems, it is probable that peatlands depend not only on the environmental conditions related to direct inputs from precipitation, snowmelt from the adjacent mountains, and groundwater infiltration and upwelling from the current year, but also from previous years. We expect that this spatiotemporal NDVI dataset will be useful for future research devoted to understanding the relationships between the regional productivity of Andean peatlands and the large-scale climate variability underlying the patterns reported in this paper.

Conclusions
In this study, we provide the first digital inventory and spatiotemporal characterization of high Andean peatland productivity for the entire Western Altiplano. Based on an exceptional cloud-free Landsat dataset, spanning 31 years, we identified 5665 peatland units and delineated their maximum long-term extent considering the years of highest NDVI productivity. We evaluated the patterns in monthly NDVI at the regional, peatland polygon, and pixel scales. The regional NDVI signal shows that high Andean peatlands have recurrent seasonal behavior with a growing season that begins during the Southern Hemisphere mid-spring, peaks during late summer, and ends around late autumn. The regional time series is relatively stable between 1986 and 2011 and is only interrupted by two small drops and recoveries in productivity between 1989 and 1993 and 1998 and 2000. However, the most productive period of the 31-year record begins in 2011 and is maintained through to the present (2017). Our pixel-based assessment of peatland productivity shows that both reductions and gains in productivity are found across latitudinal bands, but also within individual peatland units. Despite these spatially heterogeneous patterns, a larger proportion of the total area shows increasing productivity between 2011 and2017, supporting the regional greening trend over the last seven years. While beyond the scope of the present study, further research is needed to explain the observed regional and local spatiotemporal patterns in peatland productivity and to explore the relative influence of climate factors across the western Altiplano. We believe that this high-resolution database has the potential to support peatland monitoring, improve livestock management, and to better understand the roles of local and large-scale climate on peatland productivity dynamics. Finally, this inventory can provide valuable input into local water management and conservation strategies in order to protect these fragile and unique ecosystems.