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.

Keywords:

• Climate
• land cover
• land degradation
• meteorological data
• MODIS NDVI
• multivariate analysis vegetation activity

1. 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; 2014; 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 (Lanfredi et al. 2015). 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).

2. 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₂ and CH₄ (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₂ 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.

2.1. 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 (Lanfredi et al. 2015), 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.

Figure 1. Digital elevation model (DEM) of the Basilicata region provided by the Basin Authority (projected to UTM zone 33N at reference datum WGS84). The locations of the meteorological stations are reported on the DEM. Source: (http://rsdi.regione.basilicata.it/Catalogo/srv/ita/search?hl=ita#|r_basili:28FFB196-1EAA-49B4-769D-78C030B75042)

Figure 2. CORINE Land Cover (CLC2012) map of the Basilicata region (projected to UTM zone 33N at reference datum WGS84) downloaded from the Italian National Institute for Environmental Protection and Research (ISPRA) (Source: http://www.sinanet.isprambiente.it/it/sia-ispra/download-mais/). The strong fragmentation of the study area clearly appears from the map.

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 (Simoniello et al. 2015). A dry sub-humid climate characterizes the eastern part of the region where most of the vulnerable areas are located.

3. Data

3.1. 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).

Table 1. Geographic features of the sampling stations.

Name	Code	Distance to water bodies or river (m)	Distance to urban or industrial areas (m)	Aspect	Elevation (m)
Aliano	AL	266	5146	North-east	178
Bernalda	BR	79	3882	South-west	57
Craco	CR	292	1041	South	72
Ferrandina	FE	62	389	East	351
Grottole	GR	103	1092	South-east	486
Guardia Perticara	GP	181	5348	North-east	551
Matera Nord	MT	345	729	East	223
Metaponto	ME	89	1391	South-east	9
Montalbano	MJ	52	3141	North-west	151
Montescaglioso	MS	357	7335	West	32
Nova Siri	NS	123	3551	South-west	144
Pisticci	PI	93	1922	East	189
Pisticci Scalo	PS	31	1956	North-east	49
Policoro Pantano	PP	1278	1276	Flat	4
Policoro Troyli	PT	56	2862	South-west	115
San Giorgio Lucano	SG	82	3239	South-east	456
Senise	SE	62	2156	North-west	279
Stigliano	ST	52	8247	North-east	281
Tursi	TU	137	1447	South-east	58
Villa D'Agri	VA	180	1575	West	598

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.

Table 2. Available database of meteorological data. Number of analyzed matrices.

Station	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	N. of years per station
AL		X	X	X	X		X	X	X	X	X					9
BR	X		X			X	X	X	X		X	X	X	X	X	11
CR		X	X	X	X	X	X	X	X	X	X	X	X	X	X	14
FE	X	X	X	X	X	X	X	X	X	X	X		X	X	X	14
GR		X	X		X	X	X	X	X	X	X	X	X	X	X	13
GP		X	X		X		X	X	X		X	X	X	X	X	11
MT	X	X	X	X	X	X	X	X	X	X		X	X	X	X	14
ME					X	X	X	X	X	X	X	X	X	X		10
MJ		X	X	X	X	X	X	X		X			X	X	X	11
MS	X	X	X	X	X	X	X		X	X	X	X	X	X	X	14
NS	X	X	X	X	X	X	X	X	X	X	X	X	X	X	X	15
PI	X	X			X	X	X	X	X	X	X	X	X	X	X	13
PS	X	X	X	X	X	X		X	X	X	X	X	X	X	X	14
PP	X	X	X	X	X	X	X	X	X	X	X	X	X		X	14
PT	X	X			X		X	X	X	X	X	X	X	X	X	12
SG		X	X	X		X	X	X	X		X	X	X		X	11
SE	X	X						X	X	X	X	X	X	X		9
ST	X	X	X	X	X	X	X	X			X	X	X	X	X	13
TU	X	X	X	X	X	X	X	X	X	X	X	X	X	X		14
VA	X	X			X	X	X	X		X	X	X	X		X	11
No. of consideredstations per year	5	9	10	7	10	9	11	10	10	9	9	8	10	10	9

3.2. Satellite data

The assessment of vegetation conditions during the examined period (2000–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 MOD13Q1 – vegetation 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.

3.3. 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.

4. Methods

4.1. 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 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.

4.2. 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: $$\mathit{\text{NDVI}}=\frac{({\rho }_{NIR}-{\rho }_{RED})}{({\rho }_{NIR}+{\rho }_{RED})}$$

where ${\rho }_{NIR}$ and ${\rho }_{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: $$\mathit{\text{NDVI}}\_\mathit{\text{range}}{}_p={\left[\max \right(MVC \mathit{\text{NDVI}})-\text{min}(MVC \mathit{\text{NDVI}}\left)\right]}_{2000-2014}$$

Finally, the mean value $NDV{I}_{\mathit{\text{range}}}=\mathit{\text{mean}}(NDV{I}_{\mathit{\text{range}}}{}_p)$ around the meteorological stations was calculated for different plot sizes (250, 500, 1000 and 2000 m) to perform correlations with meteorological data.

5. Results

5.1. 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.

Figure 3. Mean values of NPCI calculated for each sampling station and for each meteorological variable in the period (2000–2014): NPCI_Tmin, NPCI_Tmax, NPCI_RHmin, NPCI_RHmax and NPCI_P.

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 (D’Emilio et al. 2013) 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 (PS), Craco (CR), Matera (MT), Montalbano Jonico (MJ) and Bernalda (BR). The second sub-cluster A2 is composed by the sampling stations of Tursi (TU), Stigliano (ST), San Giorgio Lucano (SG) and Policoro Pantano (PP). The first sub-cluster B1 of the cluster B is composed by the sampling stations of Metaponto (ME), Montescaglioso (MS), Nova Siri (NS), Ferrandina (FE), Grottole (GR) and Senise (SE). The second sub-cluster B2 is composed by the sampling stations of Villa d’Agri (VA), Aliano (AL), Policoro Troyli (PT), Pisticci (PI) and Guardia Perticara (GP).

Figure 4. Dendrogram of NPCI range (r-NPCI).

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.

Table 3. Centroids of the clusters.

r-NPCI	r-NPCITmin	r-NPCITmax	r-NPCIRHmin	r-NPCIRHmax
A1	0.080*	0.081	0.227*	0.174*
A2	0.198*	0.139*	0.202*	0.172*
B1	0.087*	0.048*	0.079*	0.127*
B2	0.137*	0.069	0.076*	0.108*
Mean values	0.120	0.080	0.140	0.143

Values for r-NPCI rainfalls are not reported because irrelevant.

*Statistical significant difference between centroid values and the corresponding mean values.

We have tested these differences by means of a paired t-test (level of confidence p < 5%).

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.

5.2. 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.

Table 4. NDVI range estimated over the 15-year period per different plot sizes around the meteorological stations.

Stations	250 m	500 m	1000 m	2000 m
Aliano	0.12	0.10	0.11	0.11
Bernalda	0.12	0.10	0.13	0.14
Craco	0.12	0.10	0.08	0.08
Ferrandina	0.13	0.10	0.10	0.13
Grottole	0.14	0.15	0.10	0.12
Guardia Perticara	0.39	0.36	0.27	0.22
Matera Nord	0.29	0.27	0.21	0.16
Metaponto	0.11	0.07	0.10	0.10
Montalbano	0.21	0.20	0.16	0.12
Montescaglioso	0.12	0.11	0.11	0.11
Nova Siri	0.30	0.24	0.17	0.15
Pisticci	0.09	0.09	0.08	0.07
Pisticci Scalo	0.24	0.22	0.16	0.16
Policoro Pantano	0.30	0.28	0.23	0.20
Policoro Troyli	0.19	0.17	0.15	0.11
San Giorgio Lucano	0.27	0.23	0.17	0.15
Senise	0.08	0.08	0.09	0.10
Stigliano	0.24	0.25	0.21	0.16
Tursi	0.12	0.11	0.13	0.11
Villa D'Agri	0.18	0.14	0.14	0.12

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.

5.3. 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 threshold starts to decrease dramatically and at 2000 m the two variables appear to be completely decorrelated.

Table 5. Correlation between NDVI_range and NPCI for the relative humidity.

Parameter	NPCIRHmin	NPCIRHmax
NDVI_range (250 m)	0.55	0.48
NDVI_range (500 m)	0.51	0.54
NDVI_range (1000 m)	0.45	0.47
NDVI_range (2000 m)	–	–

Only values with significance >0.5 are reported.

6. 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.

Acknowledgements

We would thank the local agency ALSIA (Agenzia Lucana di Sviluppo e di Innovazione in Agricoltura), in particular Dr E. Scalcione, for the research material made available.

Disclosure statement

No potential conflict of interest was reported by the authors.