Greening and Browning Trends on the Pacific Slope of Peru and Northern Chile

1. Introduction

Climatic conditions are considered to be the primary determinants of vegetation patterns, mediated by a complex of locally specific variables that significantly impact ecosystem services, from watershed protection and soil conservation to land use utility and resource availability [1,2]. The Köppen-Geiger climate classification model is based on the ambition to understand what determines the types of vegetation that could be found in any given region [3,4]. It is a widely used system that divides climates into five principal groups (tropical, arid, temperate, continental, and polar), each sub-divided according to rainfall and temperature data to give thirty discrete climate zones [3]. The Köppen-Geiger system produces geographical zones strongly linked to vegetation productivity and the geographical determinants of plant diversity, which makes it particularly useful in understanding and predicting regional changes in the distribution of life zones [2,5]. It is the dynamic character of climatic conditions that produce long-term changes in vegetative distribution [6,7,8]. This complex interaction significantly impacts environmental processes that determine ecosystem services and resource availability. Globally, climate change is increasingly considered to be one of the greatest threats to environmental systems and ecosystem services, with biodiversity being particularly vulnerable [9,10,11]. Atmospheric ${\mathrm{CO}}_2$ concentration, the primary limiting factor of photosynthesis, has risen from 277ppm in 1750 to 414.35ppm in 2022 [12] and is predicted to continue rising [13,14,15]. In combination with atmospheric ${\mathrm{CO}}_2$ concentration, moisture availability and temperature are also key factors for determining vegetation productivity. While spatial changes in rainfall and temperature are highly variable, the general easing of limiting factors for photosynthesis is now considered to be causing global greening [5,16].

The climatic drivers of spatio-temporal vegetation dynamics are not well understood, especially in mountainous regions [17,18]. Arid land habitats in upland regions are considered to be more sensitive to changes in climatic variables than habitats in more humid regions, with minor changes in the climatic drivers having a proportionately greater impact in such habitats. Loss of habitat driven by climate change, especially for endemic organisms in these highly specialised environments, lends increasing justification to the monitoring of environmental change in arid mountainous regions [18], however, few regional studies of vegetation dynamics in arid upland regions have been conducted [19,20].

As in many regions of the world, environmental change in the Andes is generally associated with deforestation, long-term habitat degradation and large scale infrastructural projects [21,22]. The Pacific slope of the Andes in Peru and northern Chile consists of a useful study area for examining vegetative response to climate change for four main reasons: 1. the region is both arid and mountainous with the Andes rising rapidly from the arid coastal plain to 6000m at the continental divide, only dissected by rivers that cut deep east-west valleys running to the Pacific Ocean; 2. while the region supports relatively low levels of biodiversity, it is characterised by environmental fragility and high levels of biological endemism with four discrete centres of endemicity [23,24,25]; 3. the region is characterised by sharp east-west humidity and temperate gradients, producing distinct vegetation zones; and 4. the Pacific slope is a highly vulnerable, water-scarce region with significant pressures on natural resources from several large cities in the coastal deserts including Lima, Trujillo and Arequipa in Peru to Antofagasta and Calama in northern Chile, which are totally dependent on these scarce water resources. [26,27,28]. An understanding of the spatio-temporal responses to climate change in the region is therefore vital in predicting the impacts of climate change on ecosystem services and resource availability while conserving biological endemism [11].

Of the climatic variables that determine photosynthetic productivity, temperature, moisture availability and atmospheric ${\mathrm{CO}}_2$ concentration are widely considered the most important [29,30,31]. Temperature declines by approximately 1${}^{\circ }$C with each 100m of ascent, progressively imposing a limit to photosynthetic productivity as elevation increases. Photosynthetic ${\mathrm{CO}}_2$ uptake potential is known to decline with increasing altitude as temperature falls[32]. While ${\mathrm{CO}}_2$ concentration is more evenly distributed throughout the troposphere, variation does occur below 800m asl, but only up to 300m above any given point. This variability is related primarily to plant productivity and is enhanced in regions characterised by semi-permanent atmospheric stability as in this region [33]. Precipitation and moisture availability show high regional heterogeneity across the Pacific slope of Peru [34]. This is largely determined by seasonality, altitude, aspect and atmospheric stability and is significantly impacted by the episodic El Niño phenomena, although there appears to be a nonlinear relationship between rainfall and El Niño events [35].

The plant physiologist Wladimir Köppen was instrumental in developing the widely used Köppen-Geiger climate classification system, and motivated by his interest, as a botanist, in the climatological determinants of plant distributions [3]. Recently revised and updated, this model is based on climate data which, given its connection to plant life, can be used to analyse and predict the impact of climate change on plant distributions [2,4]. Wladimir Köppen first divided vegetation zones into five broad categories, determined by combining temperature and precipitation thresholds to give Equatorial, Arid, Warm temperate, Snow and Polar climates coded A,B,C,D,and E respectively. Each of these groups is then further divided using seasonal precipitation patterns and then again according to temperature thresholds to produce thirty distinct climate zones [2,4].

The Pacific slope of Peru and northern Chile fall into four of these climate zones, listed in Table 1 and shown in Figure 1. This includes the arid, hot deserts (coded BWh), which stretch discontinuously south along the coastal deserts, from northern Peru to northern Chile, with a climate defined by monthly average temperatures exceeding 10${}^{\circ }$C and over 50% of its precipitation falling in a short wet season. Highly sensitive floristic assemblages in this zone are sparse or largely restricted to narrow riparian strips that cross the coastal plain. The arid, cold desert region, (coded BWk), lies above the hot deserts and extends from northern Peru to Chile. This region, as defined by Köppen, gradually expands to the south along the Pacific slope to become the dominant climatic zone at lower elevations in southern Peru and Chile. This zone differs from the hot deserts in having at least one month in which average temperatures lie below 0${}^{\circ }$C. Vegetation in this zone is also sparse, but locally more prolific in the fog oases or lomas. Above the cold desert lies the cold semi-arid steppe, (coded BSk), which extends along the Pacific slope in a narrow elevational strip that gradually ascends southwards into Chile and Bolivia. This zone differs from the cold deserts in that more than 50% of its precipitation falls outside a short wet season and is typified by dense thorn scrub and cacti. Ascending to the high Andes, the cold semi-arid steppe climate zone (coded BSk) gives way to the Polar Tundra (coded ET). This is a more simply defined climate with average monthly temperatures lying above freezing but not exceeding 10${}^{\circ }$C and being a region typically covered in open grasslands with relictual patches of Polylepis woodlands. When mapped, these climatic zones produce a north to south arrangement of parallel zones following the Pacific slope of Peru and northern Chile as shown in Figure 1 which indicates our study area and major Köppen-Geiger climate zones [3,4].

The use of time-series satellite data for remote sensing has developed rapidly since their inception, providing a powerful set of tools to monitor geo-spatial and temporal changes at the Earth’s surface [36]. A large number of studies looking at regional and global changes in photosynthetic productivity at the earth’s surface allow the analysis of vegetation dynamics [37]. A standard metric for measuring the health and abundance of vegetation is the normalized difference vegetation index (NDVI), which is based on the principle that photosynthetic activity in plants will absorb red light and reflect infrared light in healthy plants while the opposite holds for unhealthy plants where photosynthesis reduces or does not take place [38]. The NDVI is then derived as the normalized ratio of reflectance between red and infrared wavelengths for terrestrial surfaces. Where photosynthetic processes are increasing plant growth, an NDVI time series exhibits a positive trend known as ’greening’, or a negative trend and known as ’browning’ where plant growth is decreasing. The NDVI is strongly influenced by atmospheric aerosol loading, the angle of solar incidence, time of day and soil conditions, all of which variably distort values [39]. The Enhanced Vegetation Index (EVI), modifies the NDVI, optimising vegetation signals by correcting for these distortions using blue wavelengths to correct for atmospheric aerosols and a canopy background adjustment factor. The two vegetation indices together provide a robust means of analysing vegetation and have been widely used in studies of vegetation trends at a variety of scales across the globe.

Huete et al, 2006, using the EVI analysed trends across the Amazon basin, identifying counter-intuitive dry season greening, firing a debate about the veracity of the method and the variable impact of regional drivers of vegetative productivity [40]. Long term regional studies, driven by the need to understand the impacts of climate change on the vegetation dynamics in the circumpolar boreal forests and Arctic region, have also been widely undertaken [41,42,43]. At a more regional scale, Aide et al, 2019, identified geo-spatial greening patterns in the Andes, in which woody vegetation decreased between 1000-1500m, but increaded above that elevation [44]. Analyses of greening patterns in the Andes and European Alps have further demonstrated the value of NDVI and EVI for monitoring trends in the retreat of ice sheets [43,45]. More complex studies have analysed geospatial trends in biologically diverse and topographically complex regions from the Himalaya [17] and Andean cordilleras [46] to Peru [46]. While most studies have employed NDVI, the EVI is increasingly viewed as a more sensitive and robust measure of vegetative productivity [47].

Although studies have been undertaken into climate change impacts on vegetation dynamics in the highly biodiverse, but sparsely populated, Amazon basin, there are no such studies of the Pacific slope of the Andes, despite high levels of biological endemism in the region [40,46,48,49]. Given dramatic environmental gradients in the Andes of the Pacific slope, with highly structured climatic zonation related to both moisture availability and elevation from west to east and a strong north-south moisture gradient, it is likely that climatic drivers will have a highly differentiated geo-spatial impact on vegetation dynamics. A pattern further complicated by the combination of isolated sub-regions of biological endemism and distinctive vegetative assemblages together with complex patterns of anthropogenic pressures. This presents a complex interpretive challenge, requiring smaller scale studies across the entire Pacific slope of Peru and northern Chile to fully understand the impacts of climate change here. At present, there remains a gap in our knowledge, not only of the response to climatic drivers across these environmental gradients, but also of the habitats and ecosystem dynamics on the Pacific slope.

In this study we have carried out a trend analysis of MODIS time-series multi-spectral imagery from 2000 to 2020 as a proxy for vegetation productivity to identify geospatial greening and browning trends for the Pacific slope of Peru and northern Chile. We verify the statistical significance of remotely sensed data to produce a clear pattern in vegetation response and validate our observations with ground-truthing field sessions. We discuss the role of climatic drivers in the region by correlating our calculated vegetation response with known climatic drivers.

2. Materials and Methods

2.1. Study Area

This study is focused on the arid Pacific slope of South America between 6${}^{\circ }$S and 22${}^{\circ }$S, extending from Chiclayo in northern Peru south to Calama in northern Chile and bounded by the Pacific Ocean to the west and the continental divide to the east; in total, an area of 50,406 km${}^2$. Figure 1 shows the study area and major climate zones, as defined by the Köppen-Geiger climate classification model [3,4]. In our area of interest, greening data was principally analysed over four climate zones; the hot arid desert (BWh), cold arid desert (BWk), cold arid steppe (BSk), and polar tundra (ET). The precise area coverage of each climate zone is shown in Table 1. Although some of our satellite images have included tropical rainforest, we restrict our studies to the Pacific slope of the Andes, which does not contain rainforest climate zones in Peru. Within these climate zones, what we refer to as the ’greening strip’ is almost entirely confined to the hot and cold arid deserts and the cold arid steppe.

The abrupt rise of the Andes from sea level to 6786m in this region combines a unique blend of climatic and vegetation types and land use patterns from west to east and north to south [50,51]. The climate along the Pacific slope is strongly influenced by the Intertropical Convergence Zone (ITCZ), the episodic El Niño Southern Oscillation (ENSO) phenomena in combination with the cold Humboldt Current and orographic influence of the Andes, producing great spatial variation in both temperature and precipitation. Regional precipitation patterns are highly variable, especially in the context of the El Niño Southern Oscillation. However, precipitation on the arid pacific slope typically ranges from 50 mm/year in the northern coastal deserts to 0 mm/year in the south[52] At higher elevations up to the continental divide precipitation ranges up to 1000 mm/year [53]. The marked spatial variation in the distribution of precipitation is primarily due to steep west-east altitudinal gradients which consequently control the distribution of vegetation zones maintaining a significant environmental diversity.

The arid and hyper-arid coastal deserts of Peru and Chile are bisected by numerous broad river valleys and associated riparian vegetation with a series of coastal lomas or fog oases, on the coastal hills from Trujillo southwards. In the dramatic transition from the coast through the arid sub-tropics, vegetative associations grade from thorny cactus scrub to dry shrubby woodlands at 2500m [54]. In turn, this grades to a complex mix of relictual polylepis woodlands, puna grasslands and peatlands (wetlands) to the arctic-alpine zone above 5100m, rising locally to over 6500 masl in the Cordillera Blanca and Cordillera Occidental.

2.2. Remote Sensing Data

In this analysis, we used the EVI band of the Terra module of the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor as a proxy for land cover changes [55]. EVI was chosen in preference to NDVI and the Leaf Area Index (LAI) as it corrects for atmospheric distortion and land surface characteristics in a region where bare earth and high atmospheric aerosol loading are particularly evident [56]. The data set used is the MOD13Q1.006 Terra Vegetation Indices 16-Day Global 250m resolution provided by NASA LP DAAC at the USGS EROS Center [56]. This data set creates a composite image of the earth by integrating the best quality images from a 16-day period, with each pixel marked with an overall quality indicator. The MODIS image collections were accessed and pre-processed using the Google Earth Engine platform [57]. Changes to the EVI were investigated over a twenty-one year period from the beginning of 2000 to the end of 2020.

Our interest in climatically driven greening across the Pacific slope required that we differentiate direct anthropogenic effects, such as urban development and agricultural expansion, on the vegetation signal from broader biogeographical outcomes emergent from global or continental changes. It has been previously noted that the strong greening hot spots in the arid lands are primarily of human causation [46]. Specific EVI returns can indicate direct land cover change resulting from human activity and these were used in combination with positive field verification from satellite images to determine and confirm direct change. Direct anthropogenic land cover change can be both permanent or semi-permanent and is uniquely identifiable. For example, recently irrigated farmlands appear as patchworks of rectangular greening corresponding to fields being watered. Where irrigation has ceased, such rectangular areas appear patchy and browning as the photosynthetic capacity of the field diminishes unevenly. In the same way regularly shaped areas of extreme browning are indicative of more recent land cover transformation resulting from the construction of houses or factories and corresponding largely to the expansion of urban areas. Verifying that these extremes of land cover change corresponded to discrete greening or browning driven by human activity, thereby allowed broader identification and confirmation of regional or biogeographic change on the Pacific slope. Our field verification was undertaken in the departments of Lima and Moquegua in Peru in the Fortaleza, Pativilca, Cañete and Moquegua drainage basins from 2017 to 2022 supported with reference to Google Earth satellite imagery.

2.3. Data analysis workflow

Figure 2 details our workflow for data extraction and analysis. We separate our work into three main parts: the pixel-by-pixel analysis of the 20-year time period, the construction of spatially averaged time series, and statistical correlations with climate factors.

2.4. Pixel pre-processing

Raw EVI data were obtained from MODIS at 250m x 250m resolution [56]. In this study, we refer to `pixels’ as geographical points spaced approximately 250m by 250m apart. `Points’ refer to temporal data points in the time series of each pixel. Points are classified within the MODIS dataset as missing, good, marginal, snow/ice and clouds in the `SummaryQA’ band and labelled -1,0,1,2 and 3 respectively [58] Pixels where more than 29.3% of EVI data points were missing were not analysed because, below this threshold, a Theil-Sen regression can no longer be considered robust [59]. To ensure that the time series can be compared with an equal frequency, missing temporal data points were filled in using linear interpolation.

2.5. Seasonal decomposition

In order to extract yearly and periodic vegetation index fluctuations, the complete time series for each pixel was decomposed using additive decomposition, which resulted in three time series patterns: trend, seasonal, and residual components. We used the statsmodels.tsa.seasonal.seasonal_decompose python package [60]. The seasonal component was removed using a periodicity of 23 points per year, preserving the trend and residual components. The residual component was not removed at this stage because the package employs a self-reportedly naïve method of extracting the trend, leaving a result that is too smooth and removes all finer details from the time series [60].

2.6. EVI upper envelope calculation

Cloudy days and poor atmospheric conditions will negatively affect the EVI and produce a false reading of degraded vegetation health. To minimize this effect, we take the upper envelope of our entire time series. The methodology, based on an improved Savitzky-Golay filter, is largely adapted from Chen et. al [61].

We used linear interpolation to replace data points that were classified as snow/ice, clouds, and points that corresponded to missing EVI – this was defined as $N^0$. The completed series was then smoothed using Whittaker smoothing [62], a version of spline smoothing for discrete series with smoothing parameter $\lambda =10$. The smoothed series was defined as $N^{tr}$. $N^0$ and $N^{tr}$ were compared, and a new time series $N^1$ was defined where $N_i^1=\max (N_i^0,N_i^{tr})$. This removed the downward spikes where EVI is underestimated due to clouds. The absolute difference between $N^0$ and $N^{tr}$ is defined as $d_i=|N_i^0-N_i^{tr}|$ with the maximum difference calculated as $d_{\max }=\max \left(d_i\right)$. The fitting weights are calculated using

$$W_i=\left\{\begin{array}{cc}1,\hfill & \mathrm{when}{N}_i^0\ge N_i^{tr}\hfill \\ 1-d_i/d_{\max },\hfill & \mathrm{when}{N}_i^0<N_i^{tr}\hfill \end{array}\right..$$

(1)

A new, upper envelope time series is defined as the highest value between the original time series and a fitted time series using the Savitzky-Golay (Sav-Gol) fitting, we label this time series $N_i^{fit}$. We define this new envelope as $N_i^k=\max (N_i^0,N_i^{fit})$.

A fitting index was calculated using

$$F_k=\sum _{i=1}^n|N_i^k-N_i^0|\times W_i$$

(2)

Iteratively fitting $N^k$ and calculating $F_k$ led to the exit condition $F_{k-1}\ge F_k\le F_{k+1}$ i.e. minimum $F_k$ being reached. The expected behaviour of the iteration-fit curve is that F will decrease to a minimum and then increase, which typically takes 2 iterations. The fit corresponding to $F_k$ was taken. For the Savitzky-Golay fitting, combinations of window length and poly-order from m = 5,7,9,11,13,15,17 and d = 2,3,4,5,6 with the condition m > d, were tested. It was found that m = 9, d = 6 was the combination that most consistently produced this behaviour.

For atypical cases:

• if $F_0\le F_1$, i.e. the iteration-fit curve continually increased, the fit of the first iteration was taken,
• if after 10 iterations, the $F_{k-1}\ge F_k\le F_{k+1}$ condition was not met, i.e. iteration-fit curve continually decreased, the fit of the last iteration was taken.

2.7. Assessment of statistical significance

The final fit from the Savitzky-Golay filtering was smoothed again using a Whittaker smoothing parameter of $\lambda =2$. While $\lambda =10$ was chosen earlier to produce a very smooth time series representing the long term trend, $\lambda =2$ was chosen to reduce the amount of noise remaining from the construction of the upper EVI envelope without removing details in each time series. The Whittaker filter was chosen over Savitzky–Golay because the Savitzky–Golay filter is less reliable at the start and ends of the time series due to windowing effects and is therefore considered less robust [63].

The trend in our post-processed time series was assumed to be linear, owing to the relatively short temporal observation window, and the Mann-Kendall test [64,65] ($\tau$ rank correlation) was used to determine its statistical significance. The slope and intercept were determined using the Theil-Sen estimator, [66,67] with the intercept calculated using the Conover method [68]. This was done using the pymannkendall package [69]. The criteria $\tau >0.2$, $p>0.05$, $\mathrm{sign}\left(\tau \right)=\mathrm{sign}\left(\mathrm{slope}\right)$ were used to qualify significantly greening and browning pixels.

To calculate relative EVI change (${\Delta }_{\mathrm{EVI}}$), the trend was extrapolated to the first point of the time series – this is because seasonal_decompose truncates the beginning and end of the time series when returning trend and residual components.

$${\Delta }_{\mathrm{EVI}}=\frac{\mathrm{slope}\times (\mathbf{len}\left(\mathrm{time}\mathrm{series}\right)-1)}{\mathrm{intercept}-\mathrm{slope}\times \mathbf{len}\left(\mathrm{truncated}\mathrm{points}\right)}$$

(3)

Values of this extrapolated `initial EVI’ between 0 and $\pm 0.000001$ were set to $\pm 0.000001$ to prevent division by 0 errors. For our analysis, only pixels with EVI change between -200% and 200% were considered; changes outside this range were assumed to be unnatural or caused by a very weak initial EVI value, possibly due to the presence of water or ice.

2.8. Spatially-averaged time series

Direct anthropogenic land cover changes relating principally to irrigation schemes and urban expansion were excluded from the analysis. The Land Cover Classification System (LCCS) developed by the Food and Agriculture Organisation (FAO) identifies these areas which were spatially excluded from our analysis [70].

Only statistically significant greening pixels, as described in Section 2.7 were included in our time series analysis. We manually outlined a continuous area extending from northern Peru to northern Chile that exhibited high values of relative greening (approx. above 15%), which we define as the greening strip (GS). In order to compare the ecological dynamics of the greening strip with expected climatic behaviour, we divided the strip into the three climate zones that extend the area, namely the hot arid desert (BWh), the cold arid desert (BWk) and the cold arid steppe (BSk). We investigate the EVI time series for each climate zone both within the greening strip and over our entire area of study.

To construct our time series, we take the average EVI of all included pixels within a particular region for each point in time. The months corresponding to maxima and minima of the seasonal component were found and the `greenest’ and `brownest’ months was defined as the month before, during and after these maxima and minima respectively. The average EVI values over the greenest and brownest months were compared to the overall annual average.

2.9. Statistical correlations

Monthly precipitation values were obtained from CHIRPS [71] daily data monthly averages for each region. Sea Surface Temperature (SST) was obtained from Copernicus ERA-5 (reanalysed) [72] as a spatial average over a rectangular area between 6${}^{\circ }$ S to 30${}^{\circ }$ S and 70${}^{\circ }$ W and 80${}^{\circ }$ W. This area captures broad regions of the Pacific ocean affected by the Humboldt current. As with the EVI data we collected, seasonality was removed from these data sets prior to calculating a correlation. Monthly average ${\mathrm{CO}}_2$ data (as a trended series without seasonality) was obtained from NOAA [73]. Assuming that vegetation does not respond to climate variations instantly, we allowed for a delayed response of up to 1 month for precipitation and up to 12 months for SST and ${\mathrm{CO}}_2$. By adjusting this delay, we calculated the maximum correlation between each region’s EVI time series and the above data sets. Longer delay times were not considered for precipitation data because most of the study area consists of plants which respond to precipitation in a matter of weeks. It was assumed that soil moisture would have negligible effect due to the aridity of the area. The maximum correlations between EVI series of each pair of regions was calculated with delays of up to 12 months in both directions in time.

3. Results

We observe a complex spatially differentiated pattern of changes in EVI on the Pacific slope of Peru and northern Chile over a period from 2000-2020, with discrete zoning of greening and browning throughout the region. Figure 3 illustrates this patterning, and most notably a ’greening strip’ that extends along the pacific slope, an area characterised by arid and hyper-arid climatic conditions. The EVI in this ’greening strip’, highlighted in the inset of Figure 3, has increased by 15% during this period and extends from northern Peru to northern Chile and constitutes a hitherto undescribed novel phenomena.

There is an overall greening throughout the Andes from the foothills to the continental divide while in the coastal deserts there is a mix of greening and browning with significant areas that show no change. Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 show six different areas of the greening strip in greater detail. These areas span the entire greening strip from northern Peru to Northern Chile. The higher resolution images presented in these figures reveal distinct polygons of urban areas that appear bright red or irrigated fields that appear bright green or red (depending on whether they have recently been irrigated or have been abandoned) while the natural vegetation changes are at a lower colour intensity and its boundaries less clear cut.

Strikingly, there is a regional gap in the greening strip on the coastal pre-cordilleras in central Peru. See Figure 4. However, the ’greening’ would appear to be offset eastwards to the west facing foothills of the Cordillera Blanca rather than the Cordillera Negra. This does reflect the geography of the region, with the Cordillera Negra running north-south in parallel with the much taller Cordillera Blanca to the east. The Southern Cordillera Blanca is a small area of high peaks.

Atmospheric circulation drives the easterly flowing air masses over the Andes. In this region, the taller peaks of the Cordillera Blanca (nearly 7000m peaks) and Huayhuash force greater altitudinal uplift and vertical compression of airmasses than to the north or south. Air mass compression thus clears the Cordillera Negra, leaving it a relatively drier range with a noticeably more arid coastal strip. The greening is then recorded on the lower slopes of the Cordillera Blanca, southern Cordillera Blanca, upper Rio Fortaleza and Cordillera Huayhuash.

One of the more surprising results coming from our analysis is that the greening strip does not follow a constant altitude. Instead, the strip ascends from North to South. At the most northern part of the region, around the Lambayeque area, greening occurs from 170 to 580 masl. In the southern region, around Chiquicamata, the same greening occurs between 2850 and 4250 masl. A visualisation of the relative EVI change for transects of equal latitude is shown in Figure 9. The data represented here includes statistically significant trends within the area of interest defined as the greening strip. The minimum, maximum and mean altitude for the greening strip increases as the latitude progresses southward. We also observe that the climate zones defined in the Köppen-Geiger model similarly ascend southwards but that the ’greening’ strip does so more rapidly, crossing from lower elevation to higher elevation zones. The reason for this southward ascent remains unclear as the strip does not seem to follow isotherm or isohyet lines.

We compared the greening within the three main Köppen-Geiger zones in our region of interest (BWh, BWk, BSk) and compared the values of EVI change within and outside the area we define as the greening strip. Figure 10 shows a histogram of binned distributions of statistically significant pixel greening or browning. For each climate zone, the area outside the greening strip contains a diverse set of EVI trends, both positive and negative, which sharply drops off after a positive trend greater than 30%. For pixels within the greening strip, however, there are almost no statistically significant browning pixels and the distribution is much more positive with a dropoff around 50% relative EVI change.

The qualitative distribution of relative greening in Figure 10 shows that the vegetation within the greening strip behaves in a cohesive way across all three Köppen-Geiger zones. This is contrasted by the histograms comparing vegetation within a single Köppen-Geiger zone inside and outside of the greening strip.

Table 2 provides numerical correlations between the distributions represented in Figure 10. Here, we can clearly see that the areas defined by the greening strip in the cold arid desert (GS_BWk) and the greening strip in the cold arid steppe (GS_BSk) correlate more with each other than the areas of the same climate zones outside the greening strip. These results either indicate that the phenomenon contained within the greening strip does not follow previously defined climate zones or that the zones themselves need to be redefined in order to reflect this change in climatic conditions. The latter contradicts previous predictions by Beck et al. who claim that these climate zones will not change in the next 100 years [4]. The high correlation between the full greening strip (GS) and the sub-section that belongs to the cold arid desert (GS_BWk) exists because the greening strip is mostly located in that climate zone.

In addition to these previous observations, we see intense greening in the coastal Lomas of Peru and Chile. The Lomas are coastal hills with a highly distinctive fog-dependant vegetation, exhibiting high levels of biological endemism. These fog-adapted plants produce a highly distinctive habitat constrained by the moisture content of coastal fog or Garua. They are found in the coastal deserts from central Peru to northern Chile from 5${}^{\circ }$S to 30${}^{\circ }$S. We observe that the greening in the coastal Lomas as well as that on the Pacific slope of the Andes is not directly anthropogenic in origin. We see highly differentiated patterns of change below 1100m in the coastal deserts, with greening hot spots associated with agricultural expansion and browning hot spots associated with urbanisation and the expansion of the built environment. These represent direct changes in land cover characteristics driven by human activity.

4. Discussion

We observe a continuous strip of greening along the Pacific slope that ascends from the hot desert in the north to the temperate zone (2500-3500m) in the south with a 20-year increasing trend in EVI representing a newly described regional phenomena. From this, we can deduce that changes in arid natural vegetation have occurred but not what this consists of in terms of species assemblages and the relative abundance of each species.

Although it is possible to subdivide the 20-year study period into 10 or 5 years and study short term trends, this yields different results depending on how the period has been divided. The calculation method itself is skewed by the starting year- the result will be higher than real life if the starting year had relatively less vegetation. This is amplified because the study-area comprises of sparse vegetation by nature. It is possible that the observed 20-year trend is an anomaly on the time scale of hundreds of years but this is unverifiable using remote sensing technologies as we are limited by the availability of MODIS data. However, this result can form part of the basis of our interpretation of greening and browning trends in natural systems on the pacific slope for which long term study is required to monitor and confirm.

Exploring causality for this phenomenon is problematic, since as one moves south, the greening strip ascends, which would seem counter intuitive. That is, since plant productivity would be expected to decline as temperature declines with increasing altitude. Temperature being a primary limiting factor in photosynthesis.

The increase in global ${\mathrm{CO}}_2$ contributes to the globally observed greening phenomenon and could account for greening inside and outside of the defined strip but this cannot explain the shorter-term regionally-determined fluctuations in EVI.

We must therefore look to moisture availability, the third driver in vegetation growth. Figure 11 and Table 3 show the correlations between the 20-year trend in EVI within the greening strip and sea surface temperature (SST), precipitation, and global CO${}_2$ concentrations over the same time period.

When observing the correlations in Table 3, we focus on the relative correlations between each region rather than the actual values of each correlation. We note that there is more variation between the correlation of each climate zone outside the strip with precipitation, whereas inside the strip this more similar. Outside the strip, there is less correlation with SST and higher correlation with CO${}_2$ compared to inside the strip. It is important to note that although we correlated anomalies in EVI with precipitation in specific to each subdivided area, this does not consider the hydrology of the area which is beyond the scope of this paper.

When comparing the greening strip with land surface temperature (LST), we see complex patterns of warming and cooling. Visual confirmations with known agricultural sites have shown us that LST cooling is a response to vegetation growth rather than a driver. However, we cannot conclusively determine that LST is not a driver to vegetation growth in the greening strip since we also see regional warming highlighted by the Climate Change Institute’s reanalyzer [75].

Although the complex interplay between LST and greening or browning is beyond the scope of this research paper and will be the topic of a future publication, this preliminary observation indicates that plant feedback plays a role in the greening strip’s changing environment and possibly requires a redefinition or an updated boundary to the Köppen-Geiger Climate Zones.

5. Conclusions

We demonstrate that the greening strip is a statistically significant regional phenomenon. Although most likely driven by changes in the primary climatic variables that determine vegetative productivity it is counter intuitive for several reasons. First, temperatures drop as one ascends and so vegetative productivity would be expected to decline, all else being equal. This was shown using simulated photosynthetic ${\mathrm{CO}}_2$ uptake potential and altitude [32]. However, the greening strip behaves in the opposite direction, where the most intense relative greening rises in altitude as the strip progresses southward. Second, ${\mathrm{CO}}_2$ concentration is a more generalised global driver but for the most part evenly distributed throughout the atmosphere. There is some altitudinal variation below 300m that is associated with plant productivity and atmospheric stability, but the correlation of ground-level ${\mathrm{CO}}_2$ concentrations and altitude is no longer significant above 800m [33]. These previous results go against the clear delimitation of the greening strip. Third, although the novel phenomenon we describe is roughly located in the climate zones identified by the Köppen-Geiger climate classification system, this greening strip crosses from the hot arid desert, through the cold arid desert into the cold arid steppe. This evolution through the climate zones as the greening strip extends from north to south indicates that the drivers of this positive EVI trend are not locked to a single climate zone and also do not span the entirety of climate zones as they are currently defined. This holds especially true for the hot arid desert and the cold arid steppe.

In light of our findings, a comprehensive investigation into vegetation responses on the Pacific slope from northern Peru to northern Chile is essential to understand the geospatial impacts of global increases in CO${}_2$ concentrations and both atmospheric and marine temperatures changes. An unbalanced change in vegetation dynamics can affect the integrity of ecosystems and, through the use of remote sensing techniques, we have identified a region spanning from northern Peru to northern Chile that exhibits anomalous behaviour compared to global vegetation patterns. This study highlights the need for region-specific studies to understand the spatial dynamics of photosynthetic capacity, and vegetation response to changing climatic factors.