Hydrology and Earth System Sciences Discussions Estimation of Vegetation Cover Resilience from Satellite Time Series

Papers published in Hydrology and Earth System Sciences Discussions are under open-access review for the journal Hydrology and Earth System Sciences Abstract Resilience is a fundamental concept for understanding vegetation as a dynamic component of the climate system. It expresses the ability of ecosystems to tolerate disturbances and to recover their initial state. Recovery times are basic parameters of the vegetation's response to forcing and, therefore, are essential for describing realistic 5 vegetation within dynamical models. Healthy vegetation tends to rapidly recover from shock and to persist in growth and expansion. On the contrary, climatic and anthropic stress can reduce resilience thus favouring persistent decrease in vegetation activity. In order to characterize resilience, we analyzed the time series 1982–2003 of 8 km GIMMS AVHRR-NDVI maps of the Italian territory. Persistence probability of negative 10 and positive trends was estimated according to the vegetation cover class, altitude, and climate. Generally, mean recovery times from negative trends were shorter than those estimated for positive trends, as expected for vegetation of healthy status. Some signatures of inefficient resilience were found in high-level mountainous areas and in the Mediterranean subtropical ones. This analysis was refined by aggregating pixels ac-15 cording to phenology. This multitemporal clustering synthesized information on vegetation cover, climate, and orography rather well. The consequent persistence estimations confirmed and detailed hints obtained from the previous analyses. Under the same climatic regime, different vegetation resilience levels were found. In particular, within the Mediterranean subtropical climate, clustering was able to identify features with differ-20 ent persistence levels in areas that are liable to different levels of anthropic pressure. Moreover, it was capable of enhancing reduced vegetation resilience also in the southern areas under Warm Temperate sub-continental climate. The general consistency of the obtained results showed that, with the help of suited analysis methodologies, 8 km AVHRR-NDVI data could be useful for capturing details on vegetation cover activity at 25 local scale even in complex territories such as that of the Italian peninsula.


Introduction
Vegetation cover is known to play a key role in all of the processes linking land and atmosphere.It affects the energy, momentum, and hydrologic balance of the land surface thus rendering the land-atmosphere dynamical system closely coupled (for a review see Arora, 2002).
Ecosystems and climate influence one another on time scales ranging from seconds to millions of years (see, e.g., Sellers et al., 1995) so that interactions include short-term response of vegetation and soil to atmospheric processes as well as long-term evolution of ecosystems and soil structure (Pielke et al., 1998).
Anthropic activities further complicate the effects of these coupling mechanisms by exerting external forcing on the system (Wang et al., 2006).Land use change, land management, and human pressure can actually affect climate variables (e.g., precipitation and temperature) at both local and global scale (see, e.g., Pielke et al., 2002;Feddema, 2005).Nonlinear interactions could amplify (positive feedbacks) or attenuate (negative feedbacks) disturbances produced by human activities.By now, the fundamental importance of including vegetation cover in climatic and hydrological models as a dynamic component is widely recognized (see, e.g., Montaldo et al., 2005;Donohue et al., 2007).Especially in T. Simoniello et al.: Estimation of vegetation cover resilience the context of climate change, the assessment of vegetation responses to climatic stress and the evaluation of its ability to follow change and/or to recover from disturbances are fundamental.As an example, vegetation responses to warming, such as changes in phenology, leaf area index, and rooting depth, are able to influence soil moisture availability and evapotranspiration rates.Models that do not incorporate such direct biotic reactions can lead to erroneous predictions of ecosystem water balance (Zavaleta et al., 2003).
Some simple models (see, e.g., Brovkin et al., 1998;Zeng et al., 2002 and reference therein) have tried to describe the direct vegetation-rainfall relationship.Characteristic time scales are included in the models, but only the equilibrium relationships are discussed.In other words, time is normalized to the characteristic response scale of vegetation before running the models and only stationary asymptotic behaviors are described.
Nevertheless, models aiming at real forecasting could greatly benefit from the inclusion of real time scales.This is particularly true in a changing climate, when the hypothesis of long range stationarity does not hold.
The presence itself of characteristic time scales is an indication of resilience, which is a measure of the recovery rate of vegetation from environmental shock (Odum, 1983).The higher the resilience, the shorter the recovery time.Resilience is expected to be reduced in areas affected by prolonged stress of anthropic or climatic origin.Highly reduced resilience is in fact a potentially early sign of land degradation and desertification.On the contrary, in the presence of favorable conditions, vegetation is expected to grow and expand for long periods of time thus overstaying in a nonstationary state towards the climax condition (the stationary state of mature phases) (Odum, 1983).Therefore, estimating recovery times could provide realistic dynamical parameters for characterizing the actual current status of vegetation cover.If correctly related to structural, anthropic, and climatic features of the concerned territory, these times could provide useful hints on the complex links existing among these factors and vegetation within the hydrological and climate dynamics.
Undoubtedly, satellite observations are the main source of information on the land cover variability that can be used to improve the description of land properties.They provide not only vegetation cover maps to be directly assimilated into the models, but also a wealth of data concerning the dynamics of vegetation at different spatial and temporal scales to be used for improving the models themselves (Zhou et al., 2006;Pettorelli et al., 2005;Zeng et al., 2002).In particular, significant relationships have been reported between the satellite Normalized Difference Vegetation Index (NDVI, Tucker,1979) and some structural and functional characteristics of vegetation, such as biomass (e.g, Studer et al., 2007;Sellers, 1985), Leaf Area Index (LAI) (e.g, Maselli et al., 2004;Asrar et al., 1984), Net Primary Production (NPP) (e.g., Olofsson et al., 2007;Prince, 1991), and Absorbed Photosynthetic Active Radiation (APAR) (e.g, Wang et al., 2004;Gamon et al., 1995).
In this work, we used 8 km NDVI from the US National Oceanic and Atmospheric Administration (NOAA) -Advanced Very High Resolution Radiometer (AVHRR), for assessing characteristic recovery times of vegetation cover interannual variability in Italy for the period 1982-2003.Data were extracted from the Global Inventory Modeling and Mapping Studies (GIMMS) NDVI database.Such a dataset is widely used for characterizing vegetation at global scale or, at most, continental scale (e.g., Zhou et al., 2006;Heumann et al., 2007 and reference therein), but no application was preformed to assess its potential at higher detailed scales for locally investigating vegetation as we proposed in the current study.
The methodology we adopted was developed and previously applied to 1-km AVHRR-NDVI (Lanfredi et al., 2004).This method is devoted to the estimation of the persistence probability which, in turn, allows for estimating possible characteristic time scales (mean recovery times) related to vegetation cover resilience.In order to estimate such a persistence probability for different vegetation cover classes, we analyzed them separately according to the Corine land cover map.Orography is another source of important information for characterizing vegetation variability over the Italian territory; the Italian peninsula is mainly mountainous.Altitude differences determine both local vegetation typology and human land management practices; so, we also investigated resilience as a function of altitude.Moreover, in order to obtain a closer link between vegetation response and climate, resilience was analyzed according to local climate regions.As a final step of our investigation, we analyzed vegetation distribution on the basis of pheno-phases in order to estimate the persistence probability for the main phenological patterns.

Study area
Geographically, Italy lies in the temperate zone.It expands towards the central part of the Mediterranean Sea between 36 • N and 47 • N latitude (Fig. 1).Due to its sizable length, very irregular orography, narrow peninsular structure, and the extension of man managed areas, the Italian territory appears to be very heterogeneous both for climate and land cover variability.

Climate
The climate of the northernmost part, bordering the European continent, quite differs from that of the southernmost part, which is surrounded by the Mediterranean Sea.The Alps act as a partial barrier against north-west winds, while the Apennines protect the west side of the peninsula against winds from north-east; as a generality, the Tyrrhenian (west) coast is warmer and wetter than the Adriatic (east) coast.The great valley of the Po River has a special climate, characterized by hot summers and severe winters.In the level areas of Southern Italy and Sicily temperatures are higher.In Sardinia, conditions are more turbulent on the western side, and the island suffers from the cold winds blowing from the northwest and from the hot winds blowing from the southwest as well.

Vegetation cover
The native vegetation of Italy reflects the heterogeneity of the physical environments.At least three macro-areas of differing vegetation can be detected: the Alps, the Po Valley, and the Mediterranean-Apennine area.Native vegetation is in practice closeted in high altitude areas.In the strictly Mediterranean part of the Apennines, where the forests have been destroyed, typical scrubland vegetation, composed primarily of leathery, broad-leaved evergreen shrubs or small trees, which is called maquis, has grown up.It can be mainly found on the lower slopes of mountains bordering the Mediterranean Sea.Altogether, vegetation typologies are very heterogeneous.They range from plants typical of warm climates, such the papyrus plant that is present in Sicily, to plants that are endemic in northern Europe, such as the Norway spruce and the Scotch pine.However, the main part of the territory is man managed; agriculture is quite diffused even in non level areas.It is mainly divided into field crops, fruit tree plantations, and agro-forestry areas.Land use and land cover patterns in anthropized areas are rather erratic since urban and industrial areas coexist with cultivated and densely vegetated areas.AVHRR data at 1 km proved to be useful for characterizing such heterogeneity (Bonfiglio et al., 2002;Maselli, 2004;Simoniello et al., 2004).In this work, for the first time, the potential use of 8 km data for the same purpose has been tested.

Satellite vegetation index
The AVHRR time series analyzed in this study is the GIMMS (Global Inventory Modeling and Mapping Studies) dataset gathered from the University of Maryland.This dataset, having a spatial resolution of 8×8 km, was corrected for the principal factors (calibration instability, intercalibration, view and illumination geometry, volcanic aerosols) upsetting the AVHRR time series (Tucker at al., 2005).
Bi-weekly NDVI data from 1982 to 2003 were reprojected into the Lambert Azimuthal Equal Area Projection, sized on the study area and composed into annual MVC (Maximum Value Composite) maps (Holben, 1986) in order to obtain a synthetic parameter not influenced by seasonal patterns but suitable for estimating interannual variations of vegetation activity.The maximum of vegetation activity is one of the recognized parameters that show a strong link to climate change (Myneni et al., 1997;Heumann et al., 2007).

Ancillary data
The adopted land cover data were elaborated in the framework of the Corine Land Cover project and were collected from National Environmental Protection Agency (APAT) (http://www.clc2000.sinanet.apat.it/).The 2000 Corine map was recoded from the original 44 classes into 9 classes by taking into account the present land covers and the detectable signal at the satellite resolution.The final map (Fig. 2) contains Mixed Broad-leaved Forests, Coniferous Forests, Maquis, Sparsely Vegetated Areas and Pastures, Permanent Crops, Complex Cultivation patterns, Arable Lands, Urbanized Areas, and Water Bodies.
Analyses at different altitudes were performed by using the GTOPO30 Digital Elevation Model (Fig. 3) elaborated at USGS (http://edc.usgs.gov/products/elevation/gtopo30/gtopo30.html).The original data, derived from several sources of raster and vector topographic information, was in geographic projection with a resolution of 30 arcseconds, approximately 1 km.
Data on Italian climate regions were acquired from the Italian pedoclimatic GIS (http://www.soilmaps.it/ita/downloads.html)elaborated by the national centre for pedological cartography (CNCP).Such a recent classification (Fig. 4) was performed on long termo-pluviometric time series and identifies six main climatic regimes: Temperate sub-oceanic, Warm Temperate sub-continental, Temperate mountainous, Mediterranean sub-oceanic, Mediterranean sub-continental, Mediterranean sub-tropical, Mediterranean mountainous.All the ancillary data were reprojected in Lambert Azimuthal Equal Area and resampled at the AVHRR resolution (8 km) in order to overlap them on satellite time series.

Estimation of persistence probability
Persistence probability q(t) measures the probability that a given fluctuating field X(x, y; t ) deviates from a reference level, usually the mean value, for a time interval t before returning back to it (First-Passage probability; Redner, 2001).It enhances those collective properties that are of interest in all the problems that require the understanding of the stochastic dynamics of spatially extended objects.If we  focus on X−X r (X r as the reference level), we have to simply look at the sign of the fluctuations and q(t) is defined as the probability that such a sign has never changed since the initial time.
The technique we used in this work was based on the concept of "sign-time" distribution (Newman and Toroczkai, 1998), which, in its standard form, is substantially a histogram estimated on the surviving times of the signs of the fluctuations.
This method assumes that the given field is not noised so that all of the observed fluctuations are really generated by the dynamics we are trying to characterize.Actually, yearto-year variability of vegetation cover includes a high content of erratic variability.Thus, in order to characterize resilience on climatic scales, we focused on NDVI trends rather than on NDVI fluctuations in themselves.In practice, we studied the tendency of vegetation activity to persist in increases or decreases since the initial time so as to characterize long term behaviors.
In our adapted algorithm, the sign-time distribution is estimated from the signs of the slopes of the linear trends obtained by fitting ever longer time sequences: the first return time is defined as the first time the slope is reset to zero.This occurs when the sequence appears to be stationary and we can assume that the current status is "statistically" equivalent to the initial one (recovery).In applicative frameworks, generally we are not able to pick up the precise return time because of the discrete and noised character of observational time series; usually we detect the first passage time, when the fluctuating field has just passed the return time and the trend has changed the sign.
Let NDVI(x, y; 0) be the value of MVC-NDVI measured at first year of the considered time series in the site (x, y) and NDVI(x, y; t) be the MVC-NDVI value at the time t and site (x, y).Our algorithm (Lanfredi et al., 2004) firstly estimates linear trends over a reference initial period [0, t i ] for any pixel (x, y).The sign of any given trend, which expresses the tendency of NDVI to increase or decrease, lasts in time until inter-annual variations do not integrate destructively.Then, at the passage from one year to the next, those sites where the addition of the current inter-annual fluctuation is not able to clear the trend previously estimated or to change its sign are classified as persistent.
The surface s(x, y; t) of trend signs in a period [0, t] is defined by assigning value 1 to the pixels (x, y) where a positive trend is detected in the given period and value −1 otherwise.
At the initial time t=t i we construct the surface s(x, y; t i ).Then, we progressively add new years and compare each new surface with the previous one by keeping track of persistent sites in persistence maps P (x, y; t) that are constructed according to the rule: At the observational time T , we obtain a cumulative field: P T t i (x, y)= T t=t i+1 P (x, y; t) by adding all the persistence maps from t i+1 to T .At any pixel (x, y), the magnitude of such a field indicates the number of years during which the initial trend has not been cleared (surviving time) and its sign indicates whether the tendency was to increase or decrease.
The last step consists in estimating the function: where N (t) and N (t i ) are the number of non cleared trends having the same sign in P (x, y; t) and P (x, y; t i ), respectively.
The persistence probability is finally fitted with an exponential decay law, e − t τ , in order to estimate the mean life time (τ ).Such a type of decay has proved to be particularly effective in a previous analysis of 1 km data (Lanfredi et al., 2004).
We would underline that we did not apply any thresholding based on the trend statistical confidence.We made this choice because our aim was not to detect "change" but to study collective behaviors of vegetation activity.We tried to evaluate whether the interannual fluctuations tend to integrate constructively (positive persistence index) or destructively (negative persistence index).In the context of climatic and hydrological models, spatial correlation of vegetation resilience is an interesting major property as is time correlation.When low significance trends are aggregated in space, their spatial correlation indicates a collective behavior that cannot be considered random in any way.

Identification of Phenological Patterns
In order to group data according to their pheno-phases, we applied a clustering procedure based on the Fuzzy K-means unsupervised algorithm (see, e.g., Richard and Jia, 2006).This algorithm is based on the minimization of an objective function that evaluates the distribution of each data value with respect to the fuzzy partition.We selected such an approach since the fuzziness concept can be appropriate for detecting clusters with a homogeneous distribution.In fact, what we expect in a phenological classification is not an abrupt variation between a cluster phenology and another; instead, especially for natural vegetation, we can expect that there are gradual differences among phenological patterns linked to different altitude and climatic conditions.
From an operative point of view, once the initial cluster centroids (seeds) are chosen, the algorithm computes the degree of membership of all feature vectors in all the clusters and, then, iteratively computes new centroids until their movement is less than a predetermined threshold (1% for our analyses).Finally, a pixel is assigned to the cluster that shows the highest membership degree.Generally, such a procedure is applied on the spectral domain mainly in the context of land cover classification to group pixels having similar spectral features.In our work, we adopted such an algorithm to analyze the temporal domain, in particular, to group pixels having similar vegetation phenology, i.e. similar timing of the vegetative cycle.Phenology is a crucial component of land-atmosphere interactions (Fitzjarrald et al., 2001;Arora and Boer, 2005).The principal phenological driving factor is the type of cover, but, for the same cover, altitude and climatic conditions determine different timing in the phenophase sequence.
The identification of the main phenological patterns was performed by applying the clustering on 15-day NDVI-MVC from January to December.In order to obtain phenological responses representative of the whole time series, we computed the average values of each bi-weekly composite with a 10-year time step; in particular, the average was calculated from 1983, 1993, and 2003 data.Since the principal trainer for sub-annual NDVI dynamics is the land cover type and any deviations from the characteristic phenological pattern are mainly due to peculiar territorial conditions (e.g.climate regime and altitude), the required number of clusters as seeds for the algorithm was chosen by taking into account the number of vegetated land cover.We included also the urbanized areas since at the considered spatial resolution such pixels are made of large vegetated portions.

Results
Persistence maps were obtained for the period 1992-2003 by using 1982-1991 as a reference period.Figure 5 shows the   total persistence map P 2003 1992 (x, y) for the whole investigated area.The inspection of this map supplies some interesting information.Areas characterized by positive persistence indices are about 63% of the total.The distribution of positive and negative values seems to be rather variable over the main part of the territory, with the exception of evident extended clusters.
The main positive one is located in the upper part of the Po Valley and extends in the north-east direction.Here, most trends have very low slopes and high NDVI values.The overall situation of this area describes fluctuating vegetation activity that tends to interfere constructively across the years.On the contrary, the lower part of the Po Valley is more heterogeneous and mainly characterized by negative persistence indices.This situation describes more unstable vegetation conditions having low levels of recovery ability.The main negative clusters are located in Southern Italy and on the islands, in accordance with independent studies on areas at risk of desertification (APAT, 2006).
The persistence map reported in Fig. 5 gathers areas with different vegetation covers, orography, and local climate together.In order to estimate persistence probability and possible characteristic times, we grouped the investigated pixels according to these three features.

Characteristic times per vegetation cover class
Sample persistence probability was estimated for different vegetated Corine land cover classes and exponential best fits were computed for each class.Altogether, the mean life times reported in Table 1 describe a rather good scenario.
All the vegetated classes are characterized by positive trends that last more than the negative ones on the average.Sparsely Vegetated Areas and Pastures are an exception since they showed similar decay times (∼16 years); such a value for negative trends could be interpreted as a preliminary sign of stressed conditions for this vegetation type.The long lasting positive trends of high value vegetation (Mixed Broadleaved Forests) account for forests being able to grow and expand (Fig. 6a).Coniferous Forests (Fig. 6b) showed less efficiency but this result could be transitory since the plot of the persistence probability shows a local fast decrease around 1994 and it becomes successively slower.

Characteristic times per altitude
The analysis performed according to altitude revealed some interesting features.By sampling the altitude range (0-3743 m) at a rate of 250 m, we noted that in most cases the average recovery time from positive trends is slower than that from the negative ones.generally tends to increase to a maximum divergence of 15 years observed for the altitude band between 750 and 1000 m (Fig. 7a).The main part of the pixels belonging to this class is located on the Tyrrhenian side of the Italian peninsula and on islands.Some alpine pixels are also included.When there is a further increase in altitude, such a difference starts to decrease and successively the situation appears to be reversed until it hits the minimum (−5.8 years) that was observed for the class located at altitude >1750 m (Fig. 7b).This range covers the whole alpine zone, mainly characterized by herbaceous vegetation and a presence of some isolated Larix, P. cembra, P. mugo, and the higher part of the pre-alpine zone having coniferous forests (mainly P. sylvestris, P. cembra, L. decidua) which are widely diffused.Very few pixels belonging to this altitude range are located in the central part of the Apennine chain or in the surrounding areas of Mt.Etna in Sicily.
Figure 8, obtained by dividing altitude in classes having a range of 500 m, neatly summarizes the smooth behaviors of characteristic time scales vs. altitude.This plot is very impressive.Characteristic times of negative trends tend to slightly increase according to altitude, whereas characteristic times of positive trends decrease with altitude.The highest levels of variation are seen within positive trends, whose recovery times are more than halved during the passage from low to high altitudes.Altogether, Fig. 8   scenario, since natural vegetation covers of high ecological value are mainly confined to mountainous areas.

Characteristic times per climatic regions
Problems enhanced for high altitude seem to be confirmed by decay times estimated for the different climatic regions (Table 2).In fact, areas having a Temperate mountainous climate regime, mainly the alpine region, showed negative trends persisting for longer time periods compared to the positive ones (about 40% more); in particular, the decaying behavior of positive trends was mainly driven by the drastic trend extinction between 1997 and 1998.On the contrary, mountainous regions with a Mediterranean climate revealed the best ratio between positive and negative decay with a marked divergence between the two extinction laws (Fig. 9a).Evidence of vegetation of healthy status was found for temperate zones and for Mediterranean ecosystems under sub-oceanic climatic conditions; whereas the Mediterranean sub-continental region showed signs of stressed vegetation since the mean recovery time for negative trends is quite high (about 18 years).
As for alpine regions, sub-tropical Mediterranean areas show a reduction in photosynthetic activity that persists longer than the positive trends (Fig. 9b); such areas, located in Apulia, Sardinia, and in the southern part of Sicily, partially correspond to areas recognized as highly affected by land degradation processes (APAT, 2006).

Characteristic times per phenological patterns
Results from clustering are shown in Fig. 10, whereas Fig. 11 reports the characteristic phenology of each cluster.As can be noted, clusters related to the principal phenological patterns synthesize the main features of the land cover map, altitudes and climatic regions.For the same land cover type, the standard phenological cycle is modified (start and end timing, amplitude) by thermo-pluviometric conditions that are directly linked to altitude and climatic regime; clustering analyses are able to identify such differences.
In the northern part, a large cluster characterizes the Po Plain (cluster 4) for the most part occupied by crop cultivations.Such a cluster is also present, at times scattered, along Hydrol.Earth Syst.Sci., 12, 1053Sci., 12, -1064Sci., 12, , 2008 www.hydrol-earth-syst-sci.net/12/1053/2008/There is a net separation from the Tyrrhenian side (cluster 5) that is characterized by higher NDVI values and lower annual variation.The separation element is constituted by cluster 6 that extends along the Apennine chain and mainly represents the phytoclimatic Fagetum region covered by broadleaved forests (beeches, hornbeam, etc.) often mixed with firs.In fact, such a cluster showed the highest photosynthetic activity values for the greening period (from May to September).
The cluster corresponding to the lowest mountainous Alps and central Italy regions (cluster 7) was characterized by growing season NDVI values being generally lower than those of the forested Apennine areas, except for in the middle of the season (June-August) where they reached similar values.
The highest alpine mountains were identified by a cluster (1) showing a phenological pattern similar to the previous one, but shifted towards having lower NDVI values.
The warmest regions (Sicily, Apulia and south-eastern Basilicata) are mainly divided into two parts that generally correspond to cultivations at higher (cluster 3) and lower (cluster 2) altitudes.The hilly one shows the highest spring photosynthetic activity that rapidly decreases after the crop harvesting period; whereas the second (plains) with tree plantations presents a more constant behavior.
For each cluster, we estimated the corresponding persistence probability and the relative mean life time (Table 3).The largest difference between positive and negative trends was found for the cluster mainly corresponding to the heavy cultivated Po Plain.The negative situation highlighted within the areas under the Temperate mountainous climate was confirmed by the clustering analysis.Both the two clusters (1 and 7), identified within this climatic region, showed a recovery time from negative shocks lasting longer than that of the positive ones (Figs.12a and b).In particular, the first one is related to the alpine ecosystems at higher altitudes, whereas the second one represents the pre-alpine ecosystems with a  higher presence of broad-leaved forests and has dynamics common to the central Apennine.Several field surveys performed in these areas seem to confirm the obtained results.
In fact, a crown condition assessment, carried out for the period 1996-1999 from the CONECOFOR network (Bussotti et al., 2002), highlighted leaf color alteration and leaf damage in beech trees present in western alpine ecosystems (Piedmont) and in Norway spruces present in the middleeastern Alps (Trentino).The principal causes were mainly parasite and fungi attacks that could be indirectly linked to climate conditions; in fact, repeated dry periods accentuate the xeric condition of internal and middle Alps (Ambrosi at al., 1998).In the central Apennine region, leaf damage due to drought and fungi attacks were found for turkey oaks in 1998 and 1999 (Bussotti et al., 2002); signs of desegregation processes were also identified in successive years and a reduction of the areale extension of many sensible species was estimated for further temperature increases (Petriccione, 2007).Also in Lombardy, a classification of crown transparency showed an increase of total injury between 1995 and 1999 with a particularly increasing trend for L. deciduas and P. abies (Balestrini et al., 2002).Decay times estimated for clusters 2 and 3 confirmed results obtained for the Mediterranean sub-tropical climate areas and enhanced some differences among the Mediterranean sub-continental areas showing more stressed vegetation   covers for Apulia, Sicily and south-eastern Basilicata than for the areas in eastern Sardinia that belong to cluster 5.Such areas are in accordance with areas identified as being more sensitive to desertification at regional scale (APAT, 2006).Figure 12c shows the best fits for cluster 2 that are similar to those estimated for cluster 3.

Conclusions
Vegetation cover resilience for the Italian territory was analyzed by estimating persistence probability and characteristic decay times of NDVI trends derived from the 8 km GIMMS dataset.Values for positive and negative trends were compared in order to derive information on vegetation status, since the presence of negative photosynthetic activity trends that exhibit longer persistence than the positive ones can be interpreted as sign of reduced vegetation resilience.Analyses were performed by considering the main controlling factors of vegetation activity, i.e. type of cover, altitude, and climate.
Results obtained per cover types showed that, by averaging information on the whole peninsula, the aggregate vegetation status is quite good (all vegetated classes are characterized by positive trends that persist longer than the negative ones).
Decay times estimated per altitude ranges revealed an interesting behavior since the persistence of negative NDVI trends increases according to elevation, whereas the positive ones have an opposite relationship.This contrasting dependence on elevation suggested inefficient resilience of vegetation activity in areas at high elevation.
The analysis of mountainous regions on a climate basis showed a net differentiation between high-altitude areas under Temperate and Mediterranean regimes respectively.In particular, the Temperate mountainous areas are characterized by the longest negative trend persistence, whereas those located under the Mediterranean regime are characterized by the shortest one.This result showed that most of the vegetation located on the southern Apennines is more resilient than that belonging to alpine and pre-alpine ecosystems, which, therefore, can be considered as the main responsible for the inefficient resilience found at highest altitudes.Among the different climate regions, low vegetation resilience was also found for areas having a Mediterranean sub-tropical climate.These areas partially include regions that are identified as being at risk of desertification.
More detailed results, in term of spatial distribution and resilience, were obtained by estimating decay times according to phenological patterns derived by clustering intra-annual NDVI data.Such an analysis allowed for a better differentiation among the coldest and warmest climates.
In particular, it was able to separate alpine and pre-alpine ecosystems showing that both of them are characterized by recovery times from negative trends that are longer than the positive ones, but the first has more rapid extinction rates for positive trends.Moreover, clustering also highlighted that vegetation in the central Apennine areas having a Warm Temperate sub-continental climate shows signs of reduced resilience.These results seem to correspond with field network monitoring observations that have highlighted signs of forest deterioration that can be linked directly (drought) or indirectly (increase in parasite and fungi attacks) to climate change.
Within the warmest climates, we also found different ecosystem responses; in particular, the vegetation resilience level of the corresponding clusters are more in accordance with the level of sensitivity to desertification of these areas than the resilience levels obtained by only considering the climatic regime.This is mainly due to the efficiency of NDVI in detecting peculiar vegetation activity in anthropized areas.In fact, in such areas, the reduction in biomass production has a double component: climatic (repeated drought periods) and anthropogenic (unsuitable land management).
On the whole, the combined use of clustering based on phenological patterns and NDVI trend persistence estimations seems to be a promising tool for deriving vegetation characteristic time scales that could be used to improve hydrological and climate modeling.Moreover, the 8 km GIMMS NDVI data seem appropriate also for characterizing vegetation cover in studies at spatial scales higher than the commonly investigated global scale.

Fig. 1 .
Fig.1.Italian territory with the regional administrative boundaries and the nam locations cited in the text.

Fig. 1 .
Fig. 1.Italian territory with the regional administrative boundaries and the name of principal locations cited in the text.
period 1992-2003 obtained from GIMMS dataset.The persistence index is expressed in number of years starting from the reference period is 1982-1991.

Fig. 5 .
Fig. 5.Total persistence map of NDVI trends, P T t i (x, y), for the period 1992-2003 obtained from 8 km GIMMS dataset.The persistence index is expressed in number of years starting from 1992; the reference period is 1982-1991.

Fig. 8 .
Figure8, obtained by dividing altitude in classes having a range of 500 m, neatly summarizes the smooth behaviors of characteristic time scales vs. altitude.This plot is very impressive.Characteristic times of negative trends tend to slightly increase according to altitude, whereas characteristic times of positive trends decrease with altitude.The highest levels of variation are seen within positive trends, whose recovery times are more than halved during the passage from low to high altitudes.Altogether, Fig.8depicts a warning

Fig. 8 .
Fig. 8. Mean life times of negative and positive NDVI trends vs. altitude.

Fig. 9 .
Fig. 9. Best fits of the persistence probability q(t) of negative and positive NDVI trends obtained by exponential decay laws for climate regions: (a) Mediterranean mountains; (b) Mediterranean sub-tropical.

Fig. 9 .
Fig. 9. Best fits of the persistence probability q(t) of negative and positive NDVI trends obtained by exponential decay laws for climate regions: (a) Mediterranean mountains; (b) Mediterranean sub-tropical.

Fig. 10 .
Fig. 10.Clusters representative of the main phenological patterns obtained by applying fuzzy k-means algorithm to the intra-annual NDVI averaged from 1983-1993-2003 data.
the persistence probability q(t) of negative and positive NDVI trends obtained by ecay laws for phenological patterns: (a) cluster 1; (b) cluster 7; (c) cluster 2.

Fig. 12 .
Fig. 12. Best fits of the persistence probability q(t) of negative and positive NDVI trends obtained by exponential decay laws for phenological patterns: (a) cluster 1; (b) cluster 7; (c) cluster 2.

Table 1 .
Percentage of pixels for the reference period, mean lifetimes τ (years) and determination coefficients R 2 of the exponential best fits obtained for vegetated Corine land covers.The determination coefficient is the percentage of the sample variance that is explained by the fit.It is estimated as: R 2 = − Nt ) 2 , where N t is the count of survived trends at the year t, Nt is the mean value of the counts, and Nt is the estimates of N t supplied by the fit.

Table 2 .
Same as Table 1 but obtained for the different Italian climatic regions.

Table 3 .
Same as Table 1 but obtained for the different clusters deriving from the fuzzy k-means unsupervised classification.