Introduction
Globally more than 807 droughts were recorded between 1900 and 2004, impacting more than 1.8 billion people, with 11 million deaths, and billions of dollars in economic losses (
Below et al, 2007;
UNESCO, 2012). Due to their multiple impacts on global agricultural, hydrological, eco-environmental, and social-economical systems, droughts have been categorized the second most geographically widespread hazard after floods according to the United Nations Educational Scientific and Cultural Organization (
UNESCO, 2012). Droughts are classified among all the natural hazards as the most devastating hazard when the number of people affected by drought is taken into considerations (
Obasi, 1994;
Keshavarz et al., 2013). From the 1970s to the early 2000s, the percentage of Earth’s land area experiencing very severe droughts doubled, according to the U.S. National Center for Atmospheric Research (
NCAR, 2005).
Degradation of natural resources has become one of the major dominating issues threating ecosystem structures, functions, and services in arid and semi-arid environments. Over 40% of the world’s land surface, with more than one billion people, is considered arid and semi-arid (UNDP/UNSO, 1997;
Reynolds et al., 2003;
Bainbridge, 2012). Arid and semi-arid vegetation community structures, patterns, functions, and services are always shaped by both climate and anthropogenic drivers at different levels (
Kaplan, 2012). Rainfall input to dry lands is the main force driving both water and primary production (
Fensholt et al., 2012;
Kaplan, 2012).
The Hopi Tribe and Navajo Nation Native American reservation areas are situated in the northeastern corner of Arizona, also known as the Four Corners region, within the lower Colorado River basin. Collectively, Hopi Tribe and Navajo Nation represent over 77,699 square kilometer of land area. Our study area has only 20 active National Weather Service Cooperative Observer Program (NWS-COOP) stations collecting hydroclimatic observations, notwithstanding the numerous gaps and inconsistencies in their records. A recent assessment of the 20 active NWS-COOP stations in the area concluded that the available hydroclimatic observations suffered from gaps and quality control problems which make them difficult to use for drought monitoring (
Garfin et al., 2007).
Seasonal rainfall for the last decades has been below average leading to a decline in range conditions, an increase in soil erosion, and a reduction in surface water flows. The consecutive dry seasons in the region have major impacts on surface water availability for livestock. The livestock production in the region has been hard hit by this ongoing drought. The main water resources for livestock in the region are the stock ponds which are supplied by surface water runoff (
Crimmins et al., 2013). Livestock production plays an important role in the region’s economy, but the decade long drought had significantly affected range conditions limiting current herd sizes. For instance; the total number of cattle on Hopi lands has declined by 60% since 1994 due to decreasing forage production (Ferguson and Crimmins, 2009).
Droughts are common in the southwest United States (Gray et al., 2003;
Cook et al., 2004;
Griffin et al., 2013), but the recent drought conditions (Figs. 1(a) and 1(b)) in the region have emerged quickly and intensely, with peak drought conditions in 2002, and have lingered for well over a decade (
Breshears et al., 2005;
Crimmins et al., 2013;
Herrmann et al., 2016). Limited long-term hydro-climatological and ecological monitoring datasets have made it difficult to track and quantify the impacts through this most recent drought, creating a challenging situation for resource managers and decision makers. The numerous gaps and inconsistencies in the hydro-climatological data make the traditional drought monitoring strategies based on precipitation observations from the existing weather stations with long periods of record not ideal for this region.
Due to the complex climatic mechanisms and limited long-term temperature and precipitation monitoring sites critical for assessing the occurrence of drought frequency and associated impacts at management scales, classic drought monitoring indices, such as Palmer drought severity index (PDSI), the crop moisture index (
Palmer, 1968;
Shafer and Dezman, 1982;
Byun and Wilhite, 1999), the surface water supply index and the standardized precipitation index (SPI) only provide general regional estimates of climate variability and frequency of drought occurrence. They provide important context and can be used to diagnose and anticipate drought conditions, but cannot capture fine-scale variations at the landscape, management scale and the complex interactions between antecedent conditions, precipitation timing and frequency and trends that can emerge in vegetation communities experiencing drought stress. In addition, several studies have reported common weaknesses related to the use of these indices to monitor drought because they rely on monthly time scale and are not sufficiently accurate to capture the beginning, end, and accumulated stress of drought. (
Sivakumar et al., 2010; Mu et al., 2013).
Because of these challenges here we are proposing to use remote sensing time series to assess the spatial patterns in the impact of the recent droughts on vegetation productivity along vegetation communities and elevation gradients of the region. Remote sensing data provide fast and accurate monitoring and tracking of land surface process at different scales from local to global (
Anyamba and Tucker, 2005;
Papeş et al., 2012;
Pôças et al., 2013). Various biophysical variables derived from satellite data, such as the Normalized Difference Vegetation index (NDVI) and land surface temperature, Leaf Area Index among other parameters, have become important datasets in aiding environmental management decisions and assessments when environmental dynamics occur at different scales from short to long-term and from local to regional (
Karnieli, 2003;
White and Nemani, 2006;
Cai and Sharma, 2010).
Productivity of natural and managed agricultural lands has been widely used as an indicator of vegetation health and seasonal dynamic (
Wessels et al., 2007;
Gamon et al., 2013;
Peng et al., 2013). The timing, magnitude and spatial patterns of vegetation dynamic from intra-seasonal to decadal scales are highly influenced by both short-term climate variability and long-term climate change. This vegetation dynamic at different spatial and temporal scales can be studied using various biophysical variables derived from remote sensing data which can quantify total biomass and characterize the unique seasonal and spectral reflectance of canopy structures and functions as well as capture spatial and temporal variation (
Ryu et al., 2010;
Horion et al., 2012).
Because of the challenges associated with the scarcity of weather and vegetation monitoring stations in the study area and the desire to better understand spatially the current status and response of its vegetation cover to droughts, remote sensing data and statistical analysis techniques provide a fast and operational research tool. Our main objective is to assess and quantify drought impacts on vegetation productivity. Results will help inform future drought monitoring and impact mitigation strategies for regional resource managers and decision makers. This proposed methodology for drought monitoring and assessment, integrates land cover, climate, and topographical data with land surface remote sensing time series observations to help: 1) evaluate rainfall controls on vegetation productivity along elevation gradients, and 2) detect any trends in vegetation productivity across the study area.
Materials and methods
Materials
Study area
The Hopi Tribe and Navajo Nation reservations are situated in the northeastern corner of the state of Arizona in the lower Colorado River basin. Collectively, the Hopi Tribe and Navajo Nation Reservations cover over 77,699 square kilometers (Fig. 3(a)). The area is characterized by cold winters and very hot summers. The average annual rainfall is less than 254 mm (Fig. 2). Climatic patterns vary from south to north across the area. The seasonal amount of rainfall varies along elevation gradients but has a distinct seasonality where it is low and moderate in the early winter, increases in February and March, and then decreasing quickly in April, with a strong monsoonal pulse between July and September. May through June is a very dry period in the area (Grahame and Sisk, 2002).Vegetation cover varies throughout the area as a result of differences in rainfall and average temperature, soil types, elevations, and land management. The area vegetative cover is dominated by three major classes: 1) grass-shrub at altitudes below 1650m consisting of sparse grassland-browse types of vegetation, 2) pinyon-juniper woodland-browse species between 1650 m and 2300 m, and 3) pine forest above 2300 m. In the areas where the dominant covers are grass-shrub and pinyon-juniper woodland, sedimentary rocks control their assemblages. Pinyon-juniper woodlands are dominated by evergreen junipers and pinyon pines and represent a widespread and diverse vegetation type that occupies plateaus, foothills, rocky outcroppings at middle elevations in the semi-arid portions of the American Southwest.
NDVI time series and Ancillary data
Multi-sensor NDVI time series data were acquired from the Vegetation Index and Phenology Laboratory at the University of Arizona (vip.arizona.edu) from 1989 to 2010. This 5.6 km seamless and sensor independent NDVI record derived from AVHRR (1981‒1999) and MODIS (2000‒2014) (
Didan, 2010;
Didan et al., 2016) provides a consistent multi-sensor time series of vegetation index. The record uses a continuity algorithm to standardize the AVHRR-NDVI time series to a MODIS-NDVI equivalent to support using the time series seamlessly. Only high quality and cloud free data are used to construct this time series and any remaining gaps are filled using a simple regression model (
Didan et al., 2016). More details about the datasets can be found on the Vegetation Index and Phenology Laboratory website: (https://vip.arizona.edu/viplab_data_explorer.php). We used the 15-day NDVI composites, which further minimizes cloud and other adverse effects, to characterize the inter-annual changes of vegetation productivity over the study period.
We used climate data from the Oregon State University Parameter-elevation Regressions on Independent Slopes Model (PRISM) (
Di Luzio et al., 2008) dataset, at 4 km spatial and monthly temporal resolution. Average monthly rainfall data were generated for the period 1989‒2010 and the long term (inter-annual) average of rainfall was computed over that period. All PRISM data were resampled to 5.6 km to match the NDVI pixel size using ArcMap software. And while PRISM gridded data is based on the very limited network of weather stations in the region it remains a unique and powerful data source and tool for spatially explicit trend analysis like the one proposed here.
Topographic and elevation data are based on the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) (
Gesch et al., 2012). The region’s elevation ranges between 821 m to 3771 m (Fig. 3(a)) corresponding to echo-zones from sparse to dense vegetation cover. The pixel size of the GDEM data was resampled to 5.6 km to match the NDVI time series using the majority resampling scheme, as opposed to simple averaging, to help preserve the dominant elevation and extremes given the complex topography of the region.
Land cover types were extracted from the 2005 North American Land Cover at 250 m spatial resolution database produced by the North American Land Change Monitoring System project (
NALCMS, 2005). The land cover of the region is fairly stable due to the little natural and anthropogenic changes. According to this classification, the area has 10 diverse land cover types (Fig. 3(b)), dominated by shrubland. The land cover data were resampled to 5.6 km to match the NDVI time series using the majority resampling scheme.
Methods
The Multi-sensor NDVI time series data were used to generate the annual NDVI based productivity parameters. The resulting productivity time series data, from 1989 to 2010, were used to characterize the extent to which the ongoing droughts has influenced and shaped the spatial temporal trends across the study area. A drought assessment approach consisting of four steps is proposed (Fig. 4). The first step was to derive the annual NDVI based productivity parameters: 1) the annual cumulative NDVI, and 2) maximum annual NDVI, both of which are good proxies of vegetation response to climate drivers, especially precipitation and temperature (
Fang et al., 2001;
Wang et al., 2003). The second step was to quantify the vegetation productivity along elevation gradients by looking at the long term average of each NDVI based productivity variable. In the third step we evaluated the relationship between the long term averages NDVI based productivity parameters and the long term average rainfall. In the fourth step we determined the directional change in spatial and temporal trends in vegetation productivity across the region.
Deriving NDVI based productivity parameters
Since the aim of this work is to look at drought impacts on vegetation productivity, we elected a simple approach of generating two spatially explicit NDVI based productivity parameters using a pre-set growing season: 1) the annual cumulative NDVI (∑NDVI, a proxy of Gross/Net Primary Production, GPP/NPP), and 2) the maximum annual NDVI. These variables were computed over the period from March to November of each year in order to minimize snow impacts on NDVI (
Zhang et al., 2004;
Delbart et al., 2006;
Shi et al., 2008;
Herrmann et al., 2016). The annual cumulative NDVI, computed spatially by integrating NDVI values from March to November of each year, is a reliable proxy of biomass production and is an integrative descriptor of ecosystem functioning that has been used extensively to characterize vegetation dynamics (Huete et al., 2008;
Boschetti et al., 2013). Empirical models and biomass harvesting studies have shown a strong correlation between the integrated NDVI and biomass productivity (
Ouyang et al., 2012;
Boschetti et al., 2013). Many studies have used these variables to capture changes in vegetation productivity across the globe (
Weiss et al., 2004; Alcaraz et al., 2006). The annual maximum NDVI has been widely used to monitor the maximum level of photosynthetic activity in the canopy (
Nieto et al., 2015). It is strongly related to vegetation density and conditions. At a given location, changes in maximum NDVI results from land-use and land-cover changes, ecosystem disturbance events, or climate stressors. We used the maximum annual NDVI to look at the maximum green biomass in a given year, as well as to capture the intra-annual climatic variability impacts on the magnitude of the annual maximum NDVI.
Vegetation production – rainfall relationship
In order to understand the response of vegetation productivity to drought, we studied their correlation in each pixel using the Ordinary Least-Squares Regression (OLS) model. This technique has been used widely to study the relationship between vegetation productivity and climatic and environmental factors (
Nicholson and Farrar, 1994;
Yuan and Roy, 2007). In order to study this relationship, we averaged spatially the annual NDVI based productivity parameters and the annual rainfall from 1989 to 2010. This analysis looked at the regression parameters by examining the strength (
p<0.05) of the coefficient of determination, (
R2), and the spatial distribution of the model residuals. This allowed evaluating the model performance and to which level vegetation productivity is controlled by rainfall distribution across the area. The residual analysis of the observed NDVI based productivity and the modeled NDVI based productivity by regression model reveals the areas where the correlation is weak and/or where the changes in vegetation productivity are controlled by more than rainfall.
Long-term trends in vegetation productivity
Parametric and non-parametric approaches are standard statistical models for assessing vegetation responses to climate change and variability, phenological changes, crop status, and land cover changes (
Wright et al., 2012;
Yin et al., 2012). A robust approach to detect changes in vegetation productivity is the use of the OLS regression model (
Liu and Gong, 2012;
Wang, 2012). This statistical technique is widely used for evaluating trends in vegetation productivity (
Ma and Frank, 2006;
Liang et al., 2012). In this study, the general trend in vegetation productivity changes was determined based on the time series analysis of the NDVI based productivity parameters. In this model, time (22 years, 1989‒2010) is the independent variable; while the dependent variables were the cumulative and maximum annual NDVI parameters. The statistical model helps identify the trends in vegetation productivity over time following this equation:
where X is the independent variable (time), Y is the modeled change variable (dependent variable), a is the slope, b presents the intercept, and e is the residual error term. Using this model, we analyzed the slopes and strength of the coefficient of determination of the regressions considering only significant change (p<0.05) to characterize the inter-annual trends in vegetation productivity. The slope of the regression serves as an indicator of the direction and strength of changes over time.
Results and discussion
Long-term average vegetation productivity
Figure 5 shows the long term NDVI based maximum vegetation productivity parameter. The spatial distribution of the long-term average of ∑NDVI and maximum NDVI for the period 1989‒2010 across the study area show that vegetation productivity increases with elevations as expected. This is attributed to- and explained by the spatial patterns of rainfall along the elevation gradients. Increase in elevation is usually associated with an increase in precipitation, cooler temperature and lower evapotranspiration demand. Deviations from this general spatial relationship indicates change in conditions, such as disturbances like fires, soil degradation, change in land management practices, or other factors. The moderating effect of rainfall along the mid and high elevations (2300‒ 3771 m) is very strong on the ∑NDVI and maximum NDVI parameters (Fig. 5). High elevation areas are dominated by high biomass coniferous forest while low biomass shrubland and grassland dominate the lower elevation.
Spatial correlation between vegetation productivity and rainfall
The relationship between NDVI based productivity and rainfall shows significant correlations, with R2 values between 0.79–0.86 for SNDVI and Max-NDVI versus rainfall. Around 91.78% of the area shows a significant relationship between the long term average SNDVI and the long term average rainfall, while 8.22% of the area shows no significant relationship. These results confirm the strong dependency of vegetation productivity on rainfall in the area, particularly for shrublands and grassland.
The residuals vary between ‒32 and 90 (VI units) across the study area at the 95% confidence level (Figs. 6(a) and 6(b)). The range between [‒32 to 0] (VI unit) corresponds to overestimation relationship, while [32 to 88] indicates underestimation relationship. Values between [0 to 90] correspond to the 95% confidence level and a strong correlation. The residuals map (Fig. 6(a)) reveals the areas where the correlation is significant or non-significant, or where the vegetation productivity is controlled by more factors than rainfall alone. It is important to note that high residuals with no significant relationship between the long term average SNDVI-rainfall were found at middle and high elevations where the ponderosa pine is the most dominant vegetation type. This indicates that rainfall input is not the only driver of vegetation productivity. While not explored in this work, the lack of strong dependence on rainfall in these areas needs to be investigated in order to know which environmental drivers control vegetation productivity of this vegetation type.
Figures 6(c) and 6(d) show the relationship between the long-term average rainfall and the long-term average Maximum-NDVI across the area. The overall residuals of the regression model range between ‒0.66 and 0.39 (VI units) and show that the range between [‒0.66 to ‒0.13] indicates the overestimated relationship, while [0.13 to 0.39] corresponds to the area where the relationship is underestimated. Values range between [‒0.12 to 0.12] refers to a significant relationship (p<0.05). As expected, not all pixels showed a significant spatial relationship between these variables.
The non-significant relationship between the long term average Max-NDVI and the long term average rainfall is found at middle and high elevations where ponderosa pine is the dominant cover. The significant correlation between Max-NDVI and rainfall was observed at low elevations. This shows that rainfall has a strong effect on vegetation productivity in the areas where grasslands and shrublands dominate. While the same cannot be said about ponderosa pine cover at high elevations where the correlation was weaker indicating other factors at play.
Spatiotemporal trends in vegetation productivity
Table 1 shows the results of this per-pixel linear trend analyses. Around 60.13% of the study area shows a significant decline in ∑NDVI, while 36% shows no significant change over the period of 1989‒2010. Most of this decline in ∑NDVI corresponds to grasslands and shrublands at low elevations, which are typically more sensitive to decline in precipitation.
A surprising increase in ∑NDVI (3.87% of the area) corresponds to middle and high-elevation areas, dominated by needle leaf trees. Whereas, long-term drought conditions are the most likely main driver of these observed declines, the increase observed at higher elevation could have resulted from warming temperatures stimulating an extended growing season (
Nemani et al., 2003;
Jolly et al., 2005), and the earlier snow melt dynamics, also associated with droughts, providing additional soil moisture early in the season during the drier periods.
Herrmann et al. (2016) show a similar divergent response along the elevation gradients of the US southwest, and hypothesized that dry years associated with higher than normal average temperature regimes may be creating favorable conditions for photosynthesis at high elevation niches that are usually very cold.
The direction and strength of the trend is captured by the coefficient of determination and the slope of the regression. In Fig. 7, the areas that show significant changes are dark-red, red, and coral colors. The remaining pixels with insignificant change are colored gray.
The spatial distribution of the slope and the coefficient of determination for ∑NDVI are shown in Figs. 7(a) and 7(b). Around 64% of the pixels show significant negative changes, compared to 36% with insignificant changes. Only 3.87% of the area shows a slight increase, mainly at higher elevations. A strong decrease in vegetation productivity (from ‒7 to ‒2 vegetation index unit, VI unit) is observed between 1800 m and 2000 m, this corresponds to areas dominated by ponderosa pine cover. The correlation coefficient (R2), in Fig. 7(b), varied between 0.18 and 0.8, with high values (from 0.5 to 0.84) mainly at low elevations where grasslands and shrublands are the dominant covers. There are few forest pixels, at high elevations that show fairly high correlation coefficient (0.6 to 0.8).
Vegetation productivity of the dominant vegetation classes, grasslands, shrublands, and pinyon-juniper woodland, and to a lesser extent ponderosa pine, have been experiencing a steady decline over the period of 1989‒2010 resulting from the ongoing droughts. Some higher elevation, ponderosa dominated forests, show increase in productivity during these droughts which might be the result of favorable temperature conditions (
Herrmann et al., 2016) associated with the still adequate precipitation, and change in snow melt dynamics, all spurred by- and associated with the droughts.
While the sum NDVI parameter tends to capture the overall climatic conditions during the full growing season, and in particular late in the season, the maximum NDVI responds mostly to the strength and patterns of the early growing season inputs like precipitation and snow melt. The maximum NDVI provides insights into how vegetation might respond to episodic precipitation even during droughts and provide a better measure of early growing season inputs. Approximately 50.04% of the study area shows significant decline with 2.4% showing significant increase, and the remaining 47.56% showing insignificant changes. The significant decline in the inter-annual maximum NDVI were observed across all vegetation types and elevation gradients.
The ‒327 to ‒100 [VI unit, scale 10,000] range corresponds to the same area that showed significant decrease in vegetation productivity. These areas are mostly located in mid-elevations which are well forested. The [‒99 to 0 VI unit] range corresponds to areas covered by shrublands and grasslands. Severe decreases in the inter-annual maximum NDVI is seen in middle elevations along forested areas. The spatial pixel-based distribution of the coefficient of determination resulting from the regression model follow the general elevation gradients. The areas with a strong decrease in maximum NDVI correspond to areas where the coefficient of determination varies between 0.40 and 0.57.
This decrease in vegetation productivity along the elevation gradients can be explained by water stress conditions within the study area that resulted from these prolonged droughts (Fig. 1(b)). Lower elevation areas are generally water-limited and drought conditions lead to a rapid response and decrease in productivity. The observed increase in vegetation productivity at higher elevation while complex to accurately attribute could have resulted from drought related consequences like changes in snow melt dynamic during drier and hotter years providing water early in the growing season to allow for longer vegetation activity and photosynthetic processes, or the temperature increase associated with the drier years may have also induced more favorable conditions at high elevation where temperature usually remains below freezing for extended periods of time (
Herrmann et al., 2016). It is also possible that below canopy vegetation cover may have benefitted from the reduced snow cover during the drier and hotter years. Nevertheless, the general results confirm vegetation conditions, as estimated through remote sensing data at mid and lower elevation, have been negatively impacted by the recent protracted drought. This is contrasted with a small number of very high elevation niches that may have seen slight improvements in vegetation conditions.
Herrmann et al. (2016) addressed the possibility of this positive response resulting from artifacts related to decreasing snow cover impacting the remote sensing dataset, and concluded that while snow cover does confound the seasonal signals the findings cannot be only explained by snow dynamics. It is however significantly clear that higher elevation vegetation cover seems to be less prone to drought conditions and whether this will continue as climate changes further remains to be seen.
Conclusions
The lands of the Hopi Tribe and Navajo Nation lie in an arid to semi-arid region making it vulnerable to climate change and anthropogenic pressure (grazing, fires, and land use changes). Few specific studies have looked at or examined the vegetative cover across this region in the context of the recent intense droughts. This effort looked at vegetation productivity responses to drought in the region measured by NDVI time series over the period of 1989‒2010. We used a simple analysis approach based on multi-sensor remote sensing NDVI time series data and statistical regression analysis to characterize the changes in the region vegetation. The results provide a first order look at these changes and can assist land managers and policy makers understand the ongoing impact of droughts on vegetation in the context of the larger climate change questions.
The region is undergoing major changes to its vegetation productivity and consequently ecosystem services. The period of 1989 to 2010, characterized by increasingly intense drought conditions, is seeing significant and widespread decline in vegetation productivity. In the Navajo nation around 74% of the area showed significant decline (p<0.05) in vegetation productivity. These declines are mostly impacting the shrubland and grassland communities of mid- elevation. At higher elevations the land cover was less prone to rainfall variability and subsequently droughts, driven by the plants deeper rooting system permitting access to more soil moisture, the larger biomass providing resilience during periods of droughts, and possibly due to changes in the temperature regime, associated with droughts, spurring more favorable conditions to plants activities during the usually freezing and long winter months, as well as changes in the snowmelt dynamics at these higher elevations. Whereas drought conditions intensifying over the study period are the main driver of the observed declines, the apparent increase of plants activities at higher elevation is a more complex response.
As droughts associated with hotter temperatures potentially become more prevalent over the region, remote areas with complex topography and vegetation cover will benefit from remotely sensed drought monitoring methods providing a near-real time and operational tool for monitoring and characterizing the ongoing responses and changes.
Higher Education Press and Springer-Verlag Berlin Heidelberg