Received 5 January 2015; accepted 22 January 2015; published 27 January 2015

ABSTRACT

This paper analyzed the temporal and spatial evolution of vegetation dynamics in various land covers in the basin of the Laguna Bustillos, Region of Cuauhtémoc, Chihuahua, Mexico. We used an NDVI time series for the months of March to April (early spring). The series was constructed from Landsat TM images for the period 1986-2011. The results show an increase of NDVI for vegetated areas, especially in conifer cover, while shrub and grassland showed a positive trend but with lower statistical significance. The increase in minimum temperatures in early spring, during the study period, was the most important factor in explaining the increase of NDVI in vegetated areas. A spatially distributed analysis shows large areas without an NDVI trend, corresponding to areas with sparse vegetation cover (degraded areas). Moreover, there are also areas with a negative trend (loss of vegetation), explained by the exploitation of trees to produce firewood which is mainly carried out by the ejidos in the region. These results help to focus human and financial resources in places where the benefit will be greatest.

Keywords:

Landsat TM, NDVI, Vegetation Dynamics, Vegetated Areas, Chihuahua, Mexico

1. Introduction

Vegetation dynamics play an important role in the evaluation of environmental processes due to the close relationship between the biosphere and global environmental parameters, which involves, among other things, the concentration of atmospheric CO₂ (Zeng et al., 1999.) , the influence of vegetation in the local water cycle (Beguería et al., 2003.) , landscape structure and diversity (Olsson et al., 2000.) , and erosion and sediment transport processes (Alatorre & Beguería, 2009; Alatorre, et al., 2011) .

Besides the weather, there are regional and local factors that also affect vegetation dynamics. For example, spatial differences in such dynamics are related to topographical conditions (Florinsky & Kuryakova, 1996; Pueyo & Beguería, 2007) , soil type (Farrar et al., 1994) , human management, land use (Stohlgren et al., 1998) , and geomorphological processes such as soil erosion (Alatorre et al., 2011; Alatorre et al., 2013) .

Several studies have identified changes in vegetation dynamics on a continental (Slayback et al., 2003; Delbart et al., 2006) , local and regional ( Andreu et al., 2007 ; Martínez-Villalta et al., 2008 ; Alatorre, et al., 2011 ; Vicente- Serrano et al., 2012 ) scale in recent decades. Most of the changes have been caused by human activity, particularly deforestation (Giglio et al., 2006; Achard et al., 2002; DeFries et al., 2002) and forest fires (Riano et al., 2007) ; moreover, rural abandonment and the consequent marginalization of some regions has contributed to the process of vegetation recovery ( Sluiter & de Jong, 2007 ; Vicente-Serrano et al., 2004; Vicente-Serrano et al., 2006a, 2006b ; Lasanta & Vicente-Serrano, 2007 ; Alatorre et al., 2011 ; Lasanta & Vicente-Serrano, 2012 ). In addition, several studies have shown that current climate trends may favor an increase in vegetation activity in different ecosystems around the world ( Myneni et al., 1998 ; Kawabata et al., 2001; Lucht et al., 2002 ) or a decrease of it (Maselli, 2004; Sarris et al., 2007; Vicente-Serrano et al., 2012) . For example, various patterns of vegetation dynamics have emerged in the Mediterranean region. In general, vegetation growth is often favored by the increase of temperature in areas where water is not the limiting factor (Andreu et al., 2007; Martínez-Villalta et al., 2008) , but this temperature increase has an opposite effect in arid and semi-arid areas (Vicente-Serrano et al., 2012) .

Systematic evaluations of vegetation activity in northwestern Mexico are scarce, although there are some works based on remote sensing. For example, Salinas-Zavala et al., (2002) studied the macro-regional effect of El Niño (ENSO) on indicators such as the Normalized Difference Vegetation Index (NDVI), and provided information for understanding the interannual variability of processes such as the increase in vegetation activity at large geogra- phic scales. Later authors such as Franklin et al. (2006) , Romo (2006) , and Bravo and Castellanos (2013) used the same index to monitor Aerial Primary Production (APP) in natural areas and areas used for livestock, discrimi- nating the effect of human activities and of the natural cycles of vegetation in this area of the country.

The results of these previous studies are restricted to very specific regions and periods (Central and Coastal Sonora in periods of one year), or to the use of satellite images (AVHRR) whose spatial resolution does not allow detailed studies. This causes significant research gaps in the region, making it impossible to: i) assess the changes in vegetation cover in recent decades; ii) detect regionaltrends in vegetation biomass; iii) study the changes in foliar activity of forest regions; iv) analyze how climatic variables (temperature and/or precipitation) and spatial patterns control aridity; v) determine the anthropogenic effect of land use; and vi) study the temporal and spatial variations of vegetation dynamics in degraded areas (gullies and erosion risk areas) where vegetation is sparse. In this regard, the objectives of this study were: i) analyze the temporal evolution of vegetation activity in vegetated areaslocated in the mountains and foothills that surround the study area using a homogenized time series of Landsat TM images (1986-2011); ii) determine which climatic variables control vegetation activity and define statistically significant temporal trends; and iii) analyze the spatial distribution exhibited by temporal trends of vegetation activity as an indicator of recovery or degradation of the vegetation cover, and quantify the effects of various topographic factors on these trends.

2. Study Area

The study area is located in the Bustillos lagoon basin, between 28˚13'19'' and 28˚59'35''N, and 106˚34'39'' and 107˚10'33''W (Figure 1(A)), with a total area of 3288 km². The basin is closed irregularly by the mountains of Pedernales, San Juan, Salitrera, Chuchupate, Sierra Azul and Rebote, so the only contribution of water is from rain. The basin has an average elevation of 2000 m, and is surrounded on the north, east, west and southwest by a set of peaks averaging 2400 m, with some peaks reaching up to 2887 m (Figure 1(B)). The National Water Commission (CONAGUA, 2010) indicates that the basin weather has an annual average rainfall of 415.7 mm, with a semi-arid,

temperate climate, and minimum and maximum average temperatures of 14.6˚C and 38˚C across the year, respectively (Figure 1(C)).

The valley bottomis mainly occupied by Phaeozem soils (Figure 1(D)), characterized by a marked accumula- tion of organic matter in the upper soil, which makes them fertile soils able to support a variety of crops and pas- tures. There are also Vertisols in the region, which are characterized by alternating between swelling and shrinking of the clays. This type of soils become hard in the dry season and plastic in the wet, making tilling very difficult except in the short transition periods between the two seasons. In general, with good management, these are very productive soils. Luvisols have great potential for a large number of crops because they are moderately weathered soils with a high degree of saturation. Leptosols predominate in the mountains and foothills. These soils are characterized by low depth (less than 30 cm) and their high content of gravel. They are unsuitable for cultivation, with very limited potential for tree crops or grass. It is best to keep them under forest cover, since their high susceptibility to erosion makes it necessary to control their use.

Finally, the information available about land tenure shows that the area is dominated by private ownership (Figure 1(E)), mainly the valley, which has the best conditions for the development of agricultural activities, while ejido properties are located in the mountains and foothills, which, due to their physiographic features, do not allow for a more intensive exploitation.

3. Data and Methodology

3.1. Selection and Preparation of the Database

To obtain a map of land cover and use, we used a Landsat TM scene (spatial resolution of 30 m) taken on October 2010, a month of the year in which there is a lower frequency of cloud cover and when perennial and seasonal crops are at their maximum development, occupying 100% of the cultivated area, which prevents confusion between bare soil and crop areas, or between irrigated meadows and natural grasslands, all of them occupying a large area of the surface under study. The image was geometrically corrected using ground control points (e.g., road intersections, airport runway intersections, bends in rivers features and the like) and the algorithm developed by Palá and Pons (1996) with the software Miramon, which takes into account the topographical distortion by incorporating a DTM.

For the classification of the Landsat TM (October 2010) scene, the atmospheric effect on the electromagnetic signal was corrected using the radiative transfer model 6S (Vermote et al., 1997) , that included external atmospheric information. Subsequently, the effect introduced by the lighting conditions was corrected to com- pensate for differences caused by uneven ground, for which we used an anisotropic or non-Lambertian topo- graphic reflectance model, which provides greater robustness than Lambertian models (Riano et al., 2003) .

For the NDVI time serieswe built a homogenized series of Landsat TM images for the period 1986-2011. The database includes 18 images, corresponding to the beginning of spring (March-April). The methods of geometric correction and topographic illumination correctionapplied to each of the images in the time series were the same as were applied to the image of October 2010. However, to homogenize the images of this time series we used the ATMOSC module of the Idrisi Kilimanjaro software and the Cos (t) model, which is an improvement over the Dark Object Subtraction model (DOS) (Eastman et al., 2004) . The Cos (t) model proposed by Chavez (1988) uses dark object subtraction to correct the haze effect and includes an estimate of transmittance, which represents the absorption by atmospheric gases and Rayleigh scattering. The DOS model assumes the atmospheric transmittance is 1.0 and spectral diffuse sky irradiance is 0.0 and the path radiance due to haze is estimated by specifying the digital number (DN) of objects that should have a reflectance or brightness of each band near of zero (e.g., deep clear lakes); however, transmittance (T) is calculated with the cosine of the zenith angle of the sun (90 - solar elevation). The data required for this method are: i) date and time of image acquisition; ii) solar and satellite angles; iii) the values of gain and offsetand the mean wavelength of the band that must be corrected. The gain and offsetvalues were incorporated into the units mW∙cm⁻²∙sr⁻¹∙um⁻¹ (milliWatts per square centimeter per steradian per micron). To check the units yielded by applying a gain and offset, we selected one of the bands in the visible to near-infrared range. Then it multiplies the highest possible image value (e.g., 255 for a byte image) by the gain and then add the offset. If the value that results is in the range of 10 - 30, the units are mW∙cm⁻²∙sr⁻¹∙um⁻¹, which is correct for use. However, if the units are a factor of 10 higher (e.g., in the hundreds), the units are W∙sr⁻¹∙um⁻¹. In this case, must shift the decimal place to the left by one digit for both the gain and the offset.

This time series of Landsat TM images was used to identify vegetation dynamics, as well as their temporal and spatial patterns, in the vegetated areas located in the mountains and foothills that surround the study area. Table 1 shows the date of each of the images used for the time series.

3.2. Classification Procedure Used

Definition of Thematic Categories and Training Areas

An important objective was to define ground truth pointsto define training areas of the Landsat TM image (October 2010), with maximum spectral heterogeneity, that represent the thematic categories present in the study area. Although the main objective of this study focuses on the analysis of vegetation dynamics as a function of vegetation activity in vegetated areas located in mountains and foothills, for theclassification algorithm it was necessary to establish a priori classesthat adequately represented the variability of land cover types presentin the

entire study area. This is because the algorithm of maximum likelihoodconsiders not only the average characte- ristics of the spectral signatureof each class, but also the covariance among classes, allowing for amore precise discrimination.

Aerial orthophotos were used in the establishmentof thematic classes in the scene,with the help of 250 ground truth points, and also in the selection of trainingareas for each thematic class. The degree of discrimination between categories was determined by the spectral signature of each of the categories and a contingency matrix generated by ERDAS 8.7 Software from spectral reflectance bands. Finally, NDVI values (October 2010) (res- caled to values between 0 and 1, where 0 corresponds to values of −1 and 1 to values of (1) were incorporated as an additional band; this was necessary to distinguish more robustly the spectral signature of each of the land covers and uses.

After verifying the adequacy of the training sample, we applied the maximum likelihood method for classi- fication. To validate the resulting classification, we calculated sensitivity and specificity statistics from a confu- sion matrix (Alatorre & Beguería, 2009) , with the help of 500 independent random points which were classified by interpretation of aerial orthophotos for verification.

3.3. NDVI Time Series (1986-2011)

NDVI time series were obtained from the homogenized series of Landsat TM images (1986-2011) with the purpose of analyzing vegetation dynamics as a function of vegetation activity in vegetated areas located in the mountains and foothills. The NDVI was calculated as (Rouse et al., 1974) :

(1)

where is the reflectance in the near infrared region of the electromagnetic spectrum and is the reflectance in the red region. NDVI is a measure of the photosynthetic capacity of the plants (Ruimy et al., 1994) and of leaf resistance to water vapor transfer (Tucker & Sellers, 1986) . However, some studies have shown a strong correlation of NDVI with the fraction of photosynthetically active radiation, vegetation biomass, green cover and leaf area index (e.g. Tucker, 1979; Tucker et al., 1981; Sellers, 1985 ). Thus, high NDVI values are indicative of high vegetation activity.

To analyze the effects of climate on vegetation activity, we obtained a database from the National Meteoro- logical Service (SMN) of Mexico; specifically, we asked data from the National Bank of Climate Data, which has historical records of the national climatological network (5000 stations), in some cases from the late last century to the present. The information used was obtained from the station in Cuauhtémoc, Chihuahua, Mexico (SMN Key: 8026), containing daily data of precipitation and daily maximum/minimum temperature, with standardized data from 1942 to 2010. The time series of total precipitation and maximum/minimum average temperatures were calculated from the normalized daily series, and in this case only the precipitation daily values were summed and the temperatures were averaged for the period immediately preceding the date of each image, thus, climatic series were calculated for the following periods prior to the date of the image: 15 days, 30 days, three months (January, February and March for March images; February, March and April for April images) and 6 months (October to March and November to April, respectively).

The topographic variables were also analyzed to assess their effect on vegetation activity (Figure 2), for this, we used a DTM with a resolution of 30m from the Mexican Continuum of Elevations (CME) which was made by the National Institute of Statistics and Geography (INEGI) and made available for download in (http://www.inegi.org.mx/geo/content/datosrelieve/continental/continuoelevaciones.aspx). The slope (%) and orientation (aspect) were obtained with the DTM (Figure 2(A) and Figure 2(B)); some studies have shown the importance of these factors to explain the recovery rates of vegetation (Pueyo & Beguería, 2007) . We also derived from the DTM the vertical dissection or potential for dissection of the geocomplex (Figure 2(C)); this map illustrates the categories of the types of relief according to the morphometric classification by levels of vertical dissection (Priego et al., 2010) derived of a DTM. Vertical dissection determines several of the features of landscape structure; on the one hand: i) the distribution of some of its components (e.g., the distribution of temperature, precipitation, somehow vegetation and partially soils and other surface materials); on the other ii) determines its capacity of association as a spatio-temporal organization (Priego et al., 2010) .

3.4. Influence of the Temporal Trends of Climatic Factors on the Temporal Variation of NDVI (1986-2011)

The existence or not of statistically significant temporal trends of NDVI has been used to detect processes of increase or decrease of vegetation activity. However, vegetation activity (and hencethe NDVI) can be affected by a number of natural factors, mostly climatic, which also undergo temporal variations. Therefore, when analyzing temporal series of NDVI it is important to account for the variance explained by those factors, in order to isolate the NDVI trends attributable to vegetation dynamics. Statistical tests based on the dependent variable alone such as the Man-Kendall’s trend detection time series of the NDVI do not identify the driving factors involved. As a consequence it is difficult to determine the existence of trends apart from those driving factors. Therefore, in order to determine the existence of temporal trends in the time series of mean NDVI values for each land cover class we performed a multivariate regression analysis against time (the year of acquisition ofthe images) and a set of climatic and astronomical covariates.

As a preliminary step we undertook a correlation analysis to determine the most appropriate time span for the climatologic time series. For the early spring images (March-April) wefound that the climatological series computed for the three months prior to the images had the greatest correlation with the NDVI time series.

As the date of acquisition of the image did not coincide year to year, which could have affected the NDVI (especially in this period of year, which is very close to the start of growing period), we also included as a covariate the Julian day of the image as a covariate. To check for temporal trends of the NDVI values that were not explained by the variability of climatic factors and the date of acquisition (Julian day) of the images, we also included the year of acquisition of the images as a covariate.

We used a backward stepwise procedure based on theAkaike’s information criterion (AIC) statistic, as imple- mented inthe function stepAIC in the MASS library of the R package forstatistical analysis (Venables & Ripley, 2002) . This function preserves for analysis only those variables that significantly explain the temporal evolution of NDVI for the different types of land cover, while the variables that do not contribute to explain NDVI values are rejected. The analysis of the results is based on: i) the goodness of fit and the degree of statistical significance of the regressions; ii) the selection of the explanatory variables; and iii) the value of the beta coefficient (standardized) to classify the variables according to their relative importance in explaining the NDVI. The presence of temporal trends in observed NDVI values that cannot be attributed to climatic or astronomical (Julian day) factors corres- pond to revegetation or degradation processes, which are correlated with the covariate “year” (and sign).

3.5. The Role of Topographic Factors and of the Spatial Distribution of Temporal Trends of NDVI (1986-2011) in Vegetated Areas

The analysis described in the previous section allowed us to determine the existence of statistically significant temporal trends in observed mean NDVI values for each land cover, and the relevance of climatic and astrono- mical factors. However, no spatial discrimination was made between areas with positive or negative trends. Thus, the next step was to repeat the multivariate analysis pixel by pixel, using the same set of climatic and astronomical (Julian day) covariates. This made it possible to obtain a map of the temporal trends of NDVI that are not explained by the covariates and, thereby, to observe the areas subject to degradation (negative trend) or recovery (positive trend) processes.

Finally, with the help of the map of the spatial distribution of the temporal trends of NDVI over the study area, we analyzed the degree of control exercised by topographic variables over trends: i) slope (%); ii) slope orientation; and iii) vertical dissection, all derived from a DTM (Figure 2).

4. Results and Discussion

4.1. Selection of Categories and Training Areas

Land cover variability of the study area comprised eight classes: human settlements (urban and rural), agriculture, apple orchards, bare soil (poor vegetation cover 0% - 30%), water bodies, grassland with scattered shrubs, shrub, and coniferous forest. The training sample was used to obtain spectral signatures foreach thematic class (Figure 3). The bare soil category was characterized by high values of brightness in all bands and greater spectral variability, and due to the absence of vegetation cover, NDVI values were the lowest, characteristics common to areas of bare soil. The categories of vegetation showed a typical spectral signature, with high values of reflectance in the bands of the infrared region (TM-4 and TM-5) and a sharp decrease towards the thermal region. In general, the spectral information shows a good discrimination between vegetation units. The inclusion of NDVI values in the spectral signatures undoubtedly helped us to make a better discrimination between categories of vegetation (Figure 3).

The contingency matrix obtained by applying the classification algorithm of maximum likelihood to the training sample indicates that all categories have success rates higher than 80% (Table 2). In the case of coni- ferous forests, they were confounded with the shrub category 16% of the time, mainly because these units are located in areas of transition between the two categories. This result would indicate that the resulting classification is very consistent.

4.2. Thematic Classification of the Images

Once the spectral separability of the different thematic units of the study area was validated, we proceeded to apply the method of maximum likelihood classification to obtain the land cover and use map (Figure 4). The use of an independent set of randomly selected pixels allowed for validating the classification model, which was very good (87.64% overall accuracy) (Table 3). The biggest confusion occurred with the category of apple orchards and bodies of water, with an error of commission of 10.53% and 0.00%, and an error of omission of 22.73% and 23.80%, respectively. The category of coniferous forest was the best discriminated from the rest of the units, with an error of commission of only 2.7%. However, if the NDVI values as another spectral band to build the spectral signature of each unit would have not been included, confusion would have been higher. The confusion between apple orchards and bodies of water was higher only when spectral information of the bands was used, which caused an overestimation of the area occupied by apple orchards, large areas of forests were confused as apple orchards, even slight confusion with some agricultural areas. There were also confusion between the categories of water bodies and bare soil when only the spectral information of the bands was used, this caused by the high turbidity of the bodies of water.

The spatial distribution of the area occupied by each category is showed in the Table 3. We observed that the areas occupied by the categories of human and agricultural settlements (urban and rural) and apple orchards are located in the bottom of the basin, where soils and topography are more suitable for human activities (see section on the study area). The bare soil category is found mainly in the transition zones between the valley and the mountains, particularly in the foothills, where intense erosion processes have resulted in continuous gully formations. The spatial distribution of the vegetation categories of grassland with scattered shrubs, and shrub and coniferous forest, suggest a gradual transition from the foothills to the highest parts of the mountains that surround the Laguna Bustillos basin (Figure 4).

4.3. Influence of Climate on the Temporal Trends of NDVI (1986-2011) in Vegetated Areas

The time series of the mean NDVI values showed a clear difference between the different land covers and useslocated in the mountains and foothills that surround the study area, with the coniferous forest category showing the highest mean values, followed by the categories of shrub, grassland with scattered shrubs and, finally, the bare soil category with the lowest mean values (Figure 5 and Figure 6). This progressive transition between the mean NDVI values of each of the categories is related to the spatial distribution of the different land covers (Figure 4 and Figure 5), where the lowest values are observed in foothills (bare soil and grassland with scattered shrubs), and the highest values in the higher parts of the mountain ranges that surround the study area (coniferous forest).

Preliminary visual inspection revealed remarkable differences in NDVI trends (Figure 6) related to the different land cover and use categories. In general, it appears that the NDVI values of vegetated areas show a positive trend, which is more evident for the category of coniferous forest. In contrast, the bare soil category shows

a dominant negative trend. Accordingly, foothill areas can be considered degraded areas (gullies and erosion risk areas) due to a sparse and incipient vegetation cover.

The multivariate regression analysis showed a good fit for the observed NDVI values, which showed very good significance for all land covers and uses (Table 4). The best model fit was obtained for the vegetation categories, particularly for the shrub category, and the worst fit was obtained for the bare soil category.

Moreover, we identified climatic variables that significantly control the temporal trends of NDVI, suggesting that the climatic conditions recorded in the past 26 years are important for explaining the evolution of vegetation activity. Average minimum temperature was the most significant explanatory factor, as shown by the values of the beta coefficients (standardized) (Table 4). The effect of minimum temperature was positive in all cases, reflecting the temporal importance of a warmer climate in late winter and early spring when plants start their growth period. This hypothesis is supported by the fact that coniferous forests showed a much lower influence from minimum temperatures than the other vegetation categories, as might be expected due to the perennial nature of needle leaves. The maximum temperature also had a positive effect on NDVI trends, although less than the minimum temperature. In contrast, accumulated rainfall did not have an effect on the temporal trends of NDVI for vegetation categories, indicating that the availability of water in the study area has no control on vegetation activity in early spring. Contrary to what was observed for vegetation categories, the bare soil category showed no correlation with any of the climatic factors.

Moreover, the acquisition date of the image (Julian day) was significant for all vegetation categories, but not for bare soil, demonstrating the relevance of the phenological state of vegetation at this time of year and the impor-

Note: The variables excluded by the stepwise procedure based on AIC (Akaike’s information criterion statistic) appear with “−”. ^*Level of significance α = 0.

tance of including this variable in studies where the NDVI is used as a measure of vegetation activity during critical periods of growth (Alatorre et al., 2010) . In this case, the acquisition date of the image presented a negative correlation with the temporal trend of NDVI recorded for vegetation categories, which means that in an earlier period higher NDVI values were registered than in later period. Early spring phenology is one of the vegetation traits that show higher response to climate (Badeck et al., 2004) . These results suggest that the early spring phenology of vegetation may have high temperature sensitivity in an increasingly warmer climate.

Once the evolution of the temporal trends of NDVI was explained by climatic and astronomical factors, it was possible to determine the existence of trends in the time series of NDVI (time variable, Table 4). Positive temporal trends of NDVI were found for vegetation categories, while a negative trend was found for the bare soil category. The strongest trend corresponded to the shrub category, with an increase of 7.21%, while for bare soil it was −3.56%, for the period 1986-2011.

The results obtained from regression analysis (Table 3) allowed us to make a more comprehensive interpreta- tion of the patterns observed in the temporal trends of NDVI (Figure 6). The apparent upward trend in NDVI for vegetation categories could be explained by a similar trend in minimum and maximum temperatures. Only in the case of the bare soil category the multivariate regression analysis indicated a negative trend that is not explained by climate variables, suggesting that in these areas there is a decrease of vegetation activity as a result of degrada- tion or loss of vegetation in the last decades, possibly related to intensive erosion processes.

These results are closely related to the changes of vegetation dynamics observed in other parts of the world. In the Western Spanish Pyrenees, Vicente-Serrano et al. (2004) found a positive temporal trend in NDVI values for forests and areas with good vegetation cover, which is associated with an increase in mean annual temperature, patterns of land abandonment and natural revegetation processes (Lasanta et al., 2007; Hill et al., 2008) . In the Central Spanish Pyrenees, Alatorre et al. (2011) studied the vegetation dynamics in areas with good vegetation cover, with erosion risk and with active erosion; the results showed that the increase in temperature and changes in cloud cover during the study period were the most important factors in explaining the evolution of the mean values of NDVI (1984-2007) observed in areas with good vegetation cover. A positive temporal trend in the mean values of NDVI for vegetation categories was found also by the present study, and it was shown that the temporal trends of minimum and maximum temperatures in the three months prior to the date of acquisition of the Landsat TM images exerted an opposite influence on the NDVI.

The fact that accumulated precipitation has no significant effect on the evolution of the mean values of NDVI could be explained by the fact that the availability of water in the study area is not a limiting factor for vegetation growth, which receives an average of about 415.7 mm∙year⁻¹, especially in summer, with 18% (75 mm) of the total annual rainfall occurring in winterperiod in which evaporation is lower, so water availability is not affected in early spring. In the study area the vegetated areas are generally distributed in relatively wet areas, in the mountains and foothills that surround the study area, where the role of precipitation may be relatively less important than that of temperature in the timing of the greenness onset. Shen et al. (2014) , demonstrated in their study as earlier- season vegetation has greater temperature sensitivity of spring phenology in northern hemisphere, where the role ofprecipitation may be relatively less important than that oftemperature, especially in forest areas.

Finally, significant negative trend was identified in the present study in bare soil category.Given the high intensity of land use in the region, these trends in the NDVI are easily attributable to human activity in the foothills of the study area, which has caused a degradation of the vegetation cover in the last decades, insomuch that currently these areas are considered bare soilwith an incipient vegetation cover. Regarding the latter possibility, studies of smaller areas have found similar temporal trends in the mean values of NDVI for vegetation, but human impact has been included as one of the explanatory factors. For example, Fuller (1998) analyzed seven years of mean NDVI values for Senegal during the period 1987-1993, and observed spatial differences in the temporal trends of NDVI with respect to agricultural and pasture management practices. Pelkey et al. (2000) reported that the creation of natural protected areas in Tanzania favored the growth of vegetation cover and biomass in areas with more intensive land use. The presence of a residual trend in NDVI values after accounting for climatic influence is regarded as evidence that other factors, such as human land use, affect the vegetation cover.

4.4. Influence of Topographic Factors and of the Spatial Distribution of NDVI Temporal Trends (1986-2011) on Vegetated Areas

The negative temporal trend in NDVI for the bare soil category could not be explained by climatic factors, suggesting the presence of degradation processes in some sectors of the study area. This possibility prompted a detailed assessment of the spatial distribution of temporal trends of NDVI to determine the presence of areas of degradation.

Figure 7 shows the spatial distribution of the negative and positive trends of NDVI in vegetated areas. There are large areas with positive temporal trends, which were explained by climatic factors (see section 4.3); however, it is also possible to see areas with negative trends, which can be associated with degradation processes of the vegetation cover. Moreover, the areas with no NDVI trend (stable) are areas that, due to their state of degradation, have not shown vegetation activity in recent decades (1986-2011).

The pixel by pixel determination of the temporal trends of NDVI allowed us to assess the occurrence of recovery and degradation processes of the vegetation cover in each of the land cover and use categories (Figure 8). In general, it can be seen that positive NDVI trends were present in all the categories analyzed. These results confirm that the coniferous forest category had greater recovery, followed, in descending order, by the categories of shrub, grassland with scattered shrubs and, finally, bare soil. Moreover, there were also negative trends in all categories, but they were more evident in the bare soil category, affecting 18% of it.

Finally, using the map of the spatial distribution of temporal trends of NDVI (Figure 7), we analyzed the degree of control exercised by topographic variables on the temporal trends of NDVI (Figure 9). The frequency distri- bution of topographic slopes (%) clearly shows that degradation processes are more pronounced in slopes of less than 10˚, while recovery processes become more significant as the degree of inclination increases. As for the orientation of the slopes, degradation processes are observed in areas with little insolation during the year (Northeast and East), while recovery processes have a greater presence on the slopes with opposite orientation to those with degradation (Southwest and West). Vertical dissection, in turn, helped establish that degradation processes have greater presence in foothill areas, while recovery processes are mainly located in the rugged areas of the mountains.

These results demonstrate that degradation processes are more pronounced in more accessible areas (slopes less than 10˚), which can be explained by the exploitation of trees to produce wood that is mainly carried out by the ejidos of the region (Figure 1(E)). Moreover, the fact that degradation processes have a greater presence on slopes with little insolation while recovery processes have a greater presence on slopes with greater insolation, is in line with other studies (Lasanta et al., 2000; Vicente-Serrano et al., 2004; Lasanta & Vicente-Serrano, 2006; Pueyo & beguería, 2007) . Finally, according to vertical dissection, degradation processes occur most frequently in foothills, which are characterized by being transition zones between different land occupations and where most human activities are performed.

5. Conclusions

This work shows the usefulness of remote sensing and GIS techniques for basic and applied geo-environmental research at basin and regional scales (areas of study between 10 and 10,000 km²). The use of supervised clas- sification techniques by the maximum likelihood method from an a priori set of categories (covers) allowed us to obtain a land cover and use map of the study area with a confidence level of 84%. The correct selection of training areas and the inclusion of NDVI allowed us to locate areas for each of the categories (replicas), resulting in a maximum variability of their spectral signatures.

The homogenized series of NDVI for the months of March-April (early spring) allowed us to analyze the spatial and temporal dynamics of vegetation activity in vegetated areas (coniferous forest, shrub and grassland with scattered shrubs) and in degraded areas (bare soil) with sparse vegetation. The results were spatially coherent and NDVI patterns were clear, coinciding with the spatial distribution of vegetation cover and land use. In summary, this study demonstrated that a significant increase in vegetation activity took place in a representative moun- tainous area of the central region of Chihuahua, Mexico, over the past 26 years, in early spring, which is largely explained by the temporal variation experienced by the minimum temperature in the study period (1986-2011). Coniferous forest and shrub are the categories that showed the greatest increase in vegetation activity, while the increase of activity in grassland with scattered shrubs has been moderate. Moreover, the active erosion and extreme environmental conditions present in the bare soils which are located mainly in foothills restricted the recovery process of vegetation during this time period. Finally, it is clear that a DTM is a useful tool for a first approach to perform a morphological exploration of the spatial distribution of temporal trends of NDVI and of recovery and degradation processes of vegetation. The slope and vertical dissection were the factors that best discriminated between recovery and degradation processes of the vegetation cover, showing that in the study period the loss of vegetation occurred more frequently in areas with lower slope and in foothills.The methodology and the results obtained in this research are a powerful tool for institutions that have responsibility for implementing environmental and land use management plans for mitigating the degradation processes of the areas vegetated, principally forests, allowing prevention efforts to be concentrated in places where the benefit will be greatest.

Acknowledgements

This research was supported by the project “Identification of eroded and risk erosion areas using remote sensing in the central region of Chihuahua” (UACJ-PTC-242; PROMEP/103・5/11/4377). Programa de Mejoramiento del Profesorado (PROMEP), de la Secretaría de Educación Pública (SEP), México. Dr. Amado, J.P., was supported by a postdoctoral scholarship from The National Council for Science and Technology of Mexico (CONACYT).

References

NOTES

^*Corresponding author.