Investigating climate variability and long-term vegetation activity across heterogeneous Basilicata agroecosystems

Abstract The Basilicata region summarizes many basic features of the biogeographic complexity characterizing Mediterranean countries. The intricate geomorphology and the long history of human management generated the current landscapes, which include both high-value ecosystems and areas prone to desertification. Preserving goods and services provided by such composite land cover mosaics poses many problems due to the interference/overlap of diverse natural and anthropic factors which make the correct selection of relevant parameters and the interpretation of observational data rather difficult. Here, we study interconnections between local climate and vegetation activity by correlating parameters characterizing the interannual statistics of the NDVI (Normalized Difference Vegetation Index), derived from satellite data, with a recently devised multivariate statistical index of meteoclimatic variability. We used a 15-year sequence of remote images concerning a set of plots located around meteorological ground stations of the central-eastern part of the region to pick up spatial structures in the vegetation–climate relationships. Our analyses were able to correlate spatial heterogeneity to variations in water exchanges between vegetation and atmosphere. This study represents a first step to improve the description of relevant processes to protect natural habitats and quality agriculture, therefore combating land degradation and climate change detrimental effects.


Introduction
One of the major concerns of Smart Communities is to balance the protection of the environment with the economic and social growth (Loperte and Cosmi 2015). Green economy partially responds to this basic requirement through the integral and sustainable exploitation of bioresources by using scientific knowledge and technological innovations as the substrate to ensure the effectiveness of the green approach.
The role played by Smart Communities is crucial to design management guidelines and to verify the successful integration of economic, environmental and social goals in the long term. In order to achieve such an integration, decision makers need factual information on relevant processes affecting natural resources, correctly seen within the heterogeneity of different biogeographic contexts.
The rich natural heritage and the large availability of bioresources and biodiversity make the Basilicata region a good place to implement green strategies in a general scenario of climate change. Not only more sustainable economic strategies could be a good opportunity for the low-income rural areas that are present in the region, but a smart approach could also be the key tool to solve or, at least, to minimize some heavy problems threatening this region. Climate change, variations in land use and anthropogenic pressure have caused land degradation and desertification phenomena rather significant in the Mediterranean countries. Several studies indicate that some areas of the Basilicata region are particularly prone to land degradation (Piccarreta et al. 2006;Basso et al. 2010;Imbrenda et al. 2013;Statuto et al. 2017) stressing on the importance of forest expansion and management to improve community resilience (Kelly et al. 2015).
Elements of such complex landscapes are often organized in small patches characterized by different levels of naturality. When studying aspects of biodiversity at the agroecosystem level, attention should be paid to both habitat diversity within one agroecosystem and agroecosystem diversity at regional and even larger scales. Many conservation actions to address local problems have to be therefore integrated in large-scale plans accounting for natural/anthropic ecosystem networks more than single local ecosystems. This demand of integrated management straightforwardly leads back to the need of a holistic approach to the monitoring useful to inform decision. In particular, long-term sustainable plans for the exploitation of natural resources need to integrate the risk linked to potential disturbing factors in the right climatic and physical context. Therefore, monitoring the link between vegetation activity and climate is one of the main preconditions to implement successful sustainability strategies.
In this work, we combined satellite vegetation data, widely used for monitoring vegetation health (Bajocco et al. 2012;Tasumi et al. 2014;Eckert et al. 2015;Ramos et al. 2015;Greco et al. 2018) with a recently developed multivariate statistical index (Di Leo et al. 2015;Giorgio et al. 2017) to study how local-scale vegetation activity is linked to climate spatial heterogeneity. Our analysis focused on the eastern part of the region which is largely used for agricultural activities and also includes areas vulnerable to land degradation and desertification . Here, preserving ecosystem goods and services for enhancing the local economy requires a better understanding of local vegetation-climate processes to reinforce resilience and combat the negative interference between climate change (warming, drought and so on) and anthropic stress (water overexploitation, intensive farming and so on).

Vegetation cover and climate: the rationale behind the work
Interactions between climate and vegetation are intrinsically bidirectional (Arora 2002). Climate controls the spatial distribution of the major vegetation types on a global scale. In turn, the presence of vegetation cover and its type alter physical characteristics of the land surface like albedo, roughness, evapotranspiration (biogeophysical mechanisms) and atmospheric gas composition, for example CO 2 and CH 4 (biogeochemical effects). The existence of a strict link between land cover annual cycles (phenology) and seasonality of meteorological variables is well known. Especially temperature plays a basic role in the seasonal timing of ecosystem phenologies (Simoniello et al. 2012), which feedback on the global climate mainly through the water and carbon cycles.
One of the many challenges for sustainability is understanding how vegetation may respond to changes in climate (e.g. due to increase in greenhouse gas concentrations) and, conversely, how changes in vegetation cover (e.g. due to human management) may affect evapotranspiration, soil moisture and, consequently, the CO 2 content in a future climate.
In the long term, considering the annual phenology as the basic temporal unit, we can look at how meteorological irregularities cumulate/vary in time investigating the contextual variability in vegetation productivity on interannual scales. Due to the complexity of the mechanisms involved, estimating the correlation of vegetation parameters with single climate variables could not be the best strategy as variations in climate produce correlated changes in multiple physical variables. As an example, analyses using ecosystem models suggest that biological systems are more sensitive to changes in energy and water balance (Schimel et al. 1997) than to temperature or precipitation in themselves (Bachelet et al. 2001).
An analysis of the coherence between climate variability and vegetation productivity could therefore benefit from an integrated approach to the problem by handling climate as a multivariate system rather than as a mere collection of independent parameters and looking at the correlation between vegetation activity and indexes representative of this multivariate structure rather than to correlation with single variables. Statistical dependence, if any, could indicate indirectly what are the likely mechanisms that drive the main vegetation-climate relationships in the long term.

Study area
Basilicata (Southern Italy) is a region with strong orographic contrasts ( Figure 1). The 46.8% of its territory is mountain, 45.2% is hilly and only 8% is flat. Land cover in this region largely reflects orography. In general, natural forested areas are present at high elevation, whereas man managed covers mostly characterize level areas ( Figure 2). The central-eastern part of the region, which is one of the warmest European areas , includes bare clay hills largely affected by erosion (Summa et al. 2007;Liberti et al. 2009) and the Metaponto alluvial plain where salinization phenomena are frequent (Satriani et al. 2012;Imbrenda et al. 2018). In spite of its vulnerability, the Metaponto plain represents the main agricultural area of the region, with intensive farming.
In the hilly central areas, the daily temperature range is very high both in summer (18-22 C) and in winter (10-15 C), and drought is severe. Along the coast, the climate is typically Mediterranean, with cool wet winters and hot dry summers. The marked bi-seasonality and the large daily temperature range expose land to erosion  and are predisposing factors for vegetation stress . A dry sub-humid climate characterizes the eastern part of the region where most of the vulnerable areas are located.

Meteorological data
Meteorological data used in this study were collected from 2000 to 2014 by the Regional Agency for Development and Innovation in Agriculture (Italian acronym -ALSIA) in 20 sampling sites of the Basilicata region ( Table 1).
The investigated area disposes of a sufficient density of meteorological stations, which are more rare in impervious mountain areas.
All the ALSIA monitoring sites are equipped with instruments for continuous and automatic measurement of agro-meteorological parameters. The examined variables are as follows: minimum air daily temperature (T min in C), maximum air daily temperature (T max in C), minimum daily relative humidity (RH min in %), maximum daily relative humidity (RH max in %) and daily precipitations (P in mm). In Table 2, for each monitoring station and for each year, available data are shown. We took into account only annual series in which the percentage of data missing is lower than 5% and in which there are more than three consecutive days of missing. The remaining data missing were filled using the mean value of nearest neighbours.

Satellite data
The assessment of vegetation conditions during the examined period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) was made using the NDVI (Normalized Difference Vegetation Index) time series derived from the MODIS (Moderate Resolution Imaging Spectroradiometer) sensor on board NASA's TERRA satellite. This index is available as product MOD13Q1vegetation indices -Collection 5, which is a composite of 16 days with 250 m of spatial resolution (available from the NASA Land Processes Distributed Active Archive Center, see https://earthexplorer.usgs.gov/). Annual time series were constructed from the initial database by considering the maximum annual values, according to the usual MVC (maximum value composite) procedure, which is adopted by the remote sensing community to obtain an estimation of the annual vegetation productivity.

Ancillary data
To identify possible vegetation covers in the vicinity of each meteorological station, we used: Aerial photographs (1:10,000, with a resolution less than 1 m) for the years 2000, 2008 and 2013 of Basilicata, all freely accessible as WMS layers in GIS environment from the Regional Geoportal of Basilicata Region (http://rsdi.regione.basilicata.it/); High-resolution Google Earth images.

Multivariate index NPCI
We analyzed H bi-dimensional annual matrices (with h = 1, … , H and H ¼ 247, see Table 2). We considered N descriptors (N ¼ 7): T m ( C) average air daily temperature, T min ( C) minimum air daily temperature, T max ( C) maximum air daily temperature,  RH m (%) average daily relative humidity, RH min (%) minimum daily relative humidity, RH max (%) maximum daily relative humidity and P (mm) daily precipitations. For each matrix $\cal {M} $ h ¼ [365 sampling days Â N descriptors], principal components analysis (PCA) was applied in order to highlight the correlation structure implicitly contained in the data set. For each PCA run, the NPCI (Normal Principal Component Index) values were calculated (Di Leo et al. 2015;Giorgio et al. 2017). NPCI estimates a standardized weight for each descriptor in the correlation structure. It allows to assess quantitatively the role of each descriptor in the different years and, at the same time, the role of the different descriptors in each year.
Moreover, for each station and for each descriptor, we calculated the range of NPCI in the investigated period. High values of NPCI indicate high variations in weight of the variable in the correlation structure, and low values of NPCI indicate the presence of a stable correlation structure over the time.
For highlighting differences among the sampling stations, we applied an unsupervised clustering algorithm to the ranges of NPCI and calculated for each station group the centroid values.

NDVI and NDVI_range
We adopted NDVI as an adequate proxy for vegetation productivity profitably used in various biogeographic regions and at different spatial scales (D'Emilio et al. 2018;Liu et al. 2018;Sulla-Menashe et al. 2018). NDVI is computed as follows: where q NIR and q RED are the reflectances at the near-infrared and red wavelengths of the adopted sensor (in this case, the index is directly available as product MOD13Q1, see Section 3.2).
We initially considered different parameters characterizing the NDVI temporal statistics in each site (mean, maximum, minimum, variability range). Our analyses led to select the range of the MVC-NDVI as the final, most explicative parameter summarizing vegetation productivity variation per each pixel p: Finally, the mean value NDVI range ¼ meanðNDVI rangep Þ around the meteorological stations was calculated for different plot sizes (250, 500, 1000 and 2000 m) to perform correlations with meteorological data.

Analysis of NPCI
We analyzed 247 yearly matrices (for each year, at least 12 matrices were available and for each sampling station at least 9 years of data were available).
For each matrix, NPCI values with the corresponding range on the investigated period were calculated. In Figure 3, the plot of mean values of NPCI of meteorological parameters for each sampling station is shown.
As shown in a previous paper (Giorgio et al. 2017), a low mean value of NPCI suggests that the descriptor has a constant weight in the examined period, whereas a high mean value of NPCI suggests that the descriptor has high variations in the correlation structure. If estimated over the whole study area, the NPCI of the precipitations has the lowest mean value (NPCI P ¼ 0.021 ± 0.004); the minimum temperature and the maximum relative humidity present lower mean values of NPCI (NPCI Tmin ¼ 0.067 ± 0.024 and NPCI RHmax ¼ 0.049 ± 0.021) with respect to the maximum temperature and minimum relative humidity (NPCI Tmax ¼ 0.253 ± 0.017 and NPCI RHmin ¼ 0.212 ± 0.012).
For highlighting differences between the sampling stations, we applied an unsupervised clustering algorithm  to the ranges of NPCI (r-NPCI) estimated over the 15-year period. Centroid values were calculated for each group of stations.
In Figure 4, we show the dendrogram of r-NPCI. This figure highlights the presence of two clusters, cluster A and cluster B. Inside each of these clusters, other two sub-clusters are present. The first sub-cluster A1 of the cluster A is composed by the sampling stations of Pisticci Scalo ( As shown in Table 3, the clustering procedure puts in evidence groups with different behaviours both for r-NPCI Tmin and r-NPCI Tmax, and r-NPCI RHmin and r-NPCI RHmax . In particular, r-NPCI Tmin is very variable in the sampling stations that are included in the A2 and B2 groups and r-NPCI Tmax is variable for the A2 sampling stations. Summarizing the results of the cluster analysis: the A1 group is characterized by statistical significant variations in r-NPCI RHmin , r-NPCI RHmax and r-NPCI Tmin ; the A2 group is characterized by high variations in all four variables; the B1 group is characterized by low variations in all four variables; the B2 group is characterized by statistical significant variations in r-NPCI RHmin , r-NPCI RHmax and r-NPCI Tmin which are opposite with respect to the A1 group.
The two macro-groups A and B can be labelled according to the weight of r-NPCI RHmin and r-NPCI RHmax in the correlation structure of the meteorological variables.

Analysis of NDVI
In Table 4, we reported the values of the NDVI_range selected to synthesize the temporal statistics of the NDVI calculated in the sampling stations.
With the help of aerial photographs and high-resolution Google Earth images, we firstly evaluated the possibility of land use/land cover changes. We did not observe dramatic phenomena requiring additional investigations. Most of the changes concerned little agricultural transformations probably linked to changes in EU's Common Agricultural Policy (CAP). Anyway, in agricultural sites, the interannual variability due to crop changes (e.g. crop rotation) is rather frequent and should introduce variability in the local microclimate as well (see Section 2). It is interesting to note that, on the average, the values of NDVI_range slightly decrease with the size of the area, thereby highlighting that significant changes are truly local. When we consider larger areas, randomness in the direction of the changes (increase/decrease) per pixel leads to a smoother temporal variability. No evident correspondence with the meteorological clustering seems to appear.

Spatial correlation analysis between NPCI range and NDVI
The correlation analysis between NDVI_range and NPCI, performed for all the meteorological variables, shows that NDVI is statistically dependent on the NPCI estimated for humidity variables (Table 5). Correlation significance at the 1000 m  We have tested these differences by means of a paired t-test (level of confidence p < 5%). threshold starts to decrease dramatically and at 2000 m the two variables appear to be completely decorrelated.

Discussion and conclusions
We performed a statistical analysis of meteorological data from ground stations and satellite data concerning vegetation activity to investigate the presence of spatial coherence in their interannual variability. The analysis of a multivariate index of meteoclimatic variability allowed to pick up two macro-groups (A and B) in the set of the analyzed sites which differ for the weight of relative humidity in the correlation structure of the meteorological parameters. The analysis of the correlation with vegetation seems to confirm the basic role of relative humidity to account for spatial diversity. As humidity is related to evapotranspiration, the main source of mutual statistical dependence seems to be the water exchange between the lowest atmosphere layer and vegetated surfaces, which therefore mostly explains spatial coherence in the long-term variability of climate and vegetation across the analyzed agroecosystems. This interpretation is also supported by the decrease of such a correlation with the size of the plots analyzed, as evapotranspiration has a truly local character. These estimates should be sensitive to effects of land cover and/or climate change which alter the current scenario and therefore are able to provide a quantitative parameter to evaluate them. As a consequence, looking at the multivariate statistics of climate jointly with local vegetation parameters may be a promising strategy to improve the description of relevant processes affecting vegetation cover heterogeneity within their proper biogeographic context. This holistic characterization could be useful to support Smart Communities in the development of plans to protect natural habitats and quality agriculture, therefore combating land degradation and climate change detrimental effects.