Land surface phenology of North American mountain environments using moderate resolution imaging spectroradiometer data

Remote Sensing of Environment 115 (2011) 1220–1233

Contents lists available at ScienceDirect

Remote Sensing of Environment j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / r s e

Land surface phenology of North American mountain environments using moderate resolution imaging spectroradiometer data Allisyn Hudson Dunn a,⁎, Kirsten M. de Beurs a,b a b

Virginia Polytechnic and State University, 131 Major Williams Hall, Blacksburg, VA 24061, United States The University of Oklahoma, 100 E. Boyd Street, Normon, OK 73019, United States

a r t i c l e

i n f o

Article history: Received 3 May 2010 Received in revised form 3 December 2010 Accepted 13 January 2011 Keywords: Land surface phenology MODIS NDVI NDII Start of season

a b s t r a c t Monitoring and understanding plant phenology is becoming an increasingly important way to identify and model global changes in vegetation life cycle events. High elevation biomes cover twenty percent of the Earth's land surface and provide essential natural resources. These areas experience limited resource availability for plant growth, development, and reproduction, and are one of the first ecosystems to reflect the harmful impact of climate change. Despite this, the phenology of mountain ecosystems has historically been understudied due to the rough and variable terrain and inaccessibility of the area. In addition, although numerous studies have used synoptically sensed data to study phenological patterns at the continental and global scales, relatively few have focused on characterizing the land surface phenology in mountainous areas. Here we use the MODIS/Terra + Aqua satellite 8-day 500 m Nadir BRDF Adjusted Reflectance product to quantify the land surface phenology. We relate independent data for elevation, slope, aspect, solar radiation, and temperature as well as longitude and latitude with the derived phenology estimates. We present that satellite derived SOS can be predicted based on topographic and weather variables with a significant R²adj between 0.56 and 0.62 for the entire western mountain range. Elevation and latitude exhibit the most significant influences on the timing of SOS throughout our study area. When examined at both the local and regional scales, as well as when accounting for aspect and temperature, SOS follows closely with Hopkins' Bioclimatic Law with respect to elevation and latitude. © 2011 Elsevier Inc. All rights reserved.

1. Introduction According to the current findings of the Intergovernmental Panel on Climate Change (IPCC), the world's surface air temperature has increased an average of 0.6 °C since the 1950s and is expected to rise another 1.4 to 5.8 °C in the 21st century (Pachauri & Reisinger, 2007). The effects of this rapid increase in global temperature can be seen in both marine and terrestrial environments and at all latitudes and longitudes through sea level rise, glacier decline, increased drought, shifts in the timing of the vegetative growing seasons, increases in gross photosynthetic activity, shifts in breeding and migration dates of fish and land animals, increased flooding, and an increase in extreme weather events (Allen & Walsh, 1996; Blyth et al., 2002; Cleland et al., 2007; Edwards & Richardson, 2004; Fagre et al., 2003; Goetz et al., 2005; Meehl et al., 2000; Myneni et al., 1997; Post, 2003; Goetz & Prince, 1996). Phenology is the study of periodic lifecycle events in plants and animals including bud burst, leaf out, reproduction, and migration, as related to the environment (Cleland et al., 2007; Lieth, 1974). Phenology ⁎ Corresponding author. Tel.: +1 540 761 7012. E-mail addresses: [email protected] (A.D. Hudson Dunn), [email protected] (K.M. de Beurs). 0034-4257/$ – see front matter © 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2011.01.005

plays an essential role in regulating the abundance and diversity of organisms and their ecological functions. Since temperate vegetation is finely tuned to the seasonality of its surroundings and responds to cumulative effects of daily weather as it fluctuates over time, its developmental stages are effective indicators of a changing climate (Cleland et al., 2007; Schwartz et al., 2006). Previous studies have shown that rapid increases in temperature and changes in precipitation brought on by current shifts in climate patterns can be linked to alterations in the timing of phenological processes including, a significant effect on the date of spring green-up, maximum photosynthetic activity, and length of the growing season, which are all critical points in plant development (Inouye, 2008; Körner, 2005; Myneni et al., 1997). If the timing of environmental cues continues to change as a result of changing weather patterns, the ability of plant and animal species to acclimatize could be in question (Cleland et al., 2007; Inouye & Wielgolaski, 2003). This strong correspondence between climate and phenological processes can be found in almost all ecosystems exhibiting seasonal weather variability, thus it is important to characterize the phenology of terrestrial environments world-wide (Parmesan & Yohe, 2003; Walker et al., 1995). However, there are few environments that are currently more affected by this changing climate than high latitude and high altitude regions. It has been shown that these boreal

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

environments are more sensitive to an increasingly warming annual temperature with individual species as well as entire ecosystems responding at different rates and varying intensities (Myneni et al., 1997; Goetz et al., 2005; Cleland et al., 2007; Goetz & Prince, 1996). Mountain and alpine ecosystems cover more than twenty percent of the Earth's land surface spanning an area from the equator to just near the poles (Inouye, 2008; Inouye & Wielgolaski, 2003; Price, 1981, 1990; Sphen & Korner, 2002; Weiss & Walsh, 2009). In addition to the commonly known characteristics of a marked topographic variation resulting in steep slopes and varied aspects, mountains are highly diverse systems in flora, fauna, and human ethnicity, and are found, at varying altitudes, on every continent. These regions provide 80% of fresh water used in today's society (Bernbaum, 1997; Blyth et al., 2002; Ives et al., 1997; Food and Agricultural Organization, 2000). High elevation environments also experience unique climate patterns aiding in the creation of niche vegetation zones which are limited in resources for plant growth, development, and reproduction (White et al., 2005). Despite the importance of mountain environments, their phenology has historically been understudied due to the rough and variable terrain and inaccessibility of the area (Bacher & Jeanneret, 1994; Ives et al., 1997; Tommervik, 2004). In recent years monitoring of remote environments has been aided by the use of satellite sensors. Imagery collected by these sensors have been important in increasing our understanding of global phenological patterns and provide one of the best forms of documenting and studying temporally continuous observations at numerous scales and spatial locations. Vegetation indices (VI) produced from reflectance data collected by the Moderate Resolution Spectroradiometers (MODIS) on Terra and Aqua or the Advanced Very High Resolution Radiometer (AVHRR) sensors have been related to the biomass of green leaf vegetation (Myneni et al., 1997; Tucker, 1979; Tucker et al., 1985; White et al., 1997) around the globe and aid the study of land surface phenology. Land surface phenology is the spatio-temporal study of the vegetated land surface as observed by synoptic sensors (de Beurs & Henebry, 2004; Morisette et al., 2009; White et al., 1997). The normalized difference vegetation index (NDVI), calculated using the near infrared (MODIS: 841 to 876 nm; AVHRR: 730 to 980 nm) and red (MODIS: 620 to 670 nm; AVHRR: 585 to 680 nm) reflectance bands, is frequently used to monitor vegetation growth cycles and health (Cleland et al., 2007; Myneni et al., 1997; Tucker, 1979). Although NDVI is the most commonly used index, it saturates at high leaf area index values and is sensitive to non-vegetation external factors, such as snow fall and snow melt that are prevalent in mountain environments. To help discriminate the effects of vegetation phenology from unrelated background phenomena such as snow, researchers have begun to expand their attention towards different vegetation indices (de Beurs & Henebry, 2010a,b; de Beurs & Townsend, 2008; Delbart et al., 2005; Delbart et al., 2006; Huete et al., 2002, 2006; Jiang et al., 2008). With the existence of the blue band, (MODIS: 459 to 479 nm; SPOT-VGT 10: 429 to 469 nm), the SWIR band (MODIS: 1230 to 1250; SPOT-VGT: 1579 to 1749 nm), and the NIR band (Spot-VGT10: 780 to 890 nm) both MODIS and SPOT-VGT S10 sensors have been used for this purpose (de Beurs & Henebry, 2010a,b; de Beurs & Townsend, 2008; Delbart et al., 2005, 2006; Huete et al., 2002; Jiang et al., 2008). Delbart et al. (2005) argue that the use of NDVI in boreal regions suffers from uncertainty in predicting start of season (SOS) due to confounding effects of snow melt. They further suggest the use of the normalized difference water index (NDWI), which is based on the near infrared and short wave infrared bands, to reduce these uncertainties. This approach had earlier been proposed by others in 2002 (Xiao et al., 2002; Xiao et al., 2002). In mountain environments where snow melt may trigger an increase in NDVI despite the lack of actual vegetation, alternative indices such as the NDWI or the normalized difference infrared index (NDII) prove useful for teasing out snowmelt versus vegetation growth. Here we compare start of season metrics based on NDVI and NDII for characterizing the land surface phenology of mountainous regions.

1221

Due to the rough terrain, harsh climate and remoteness of many mountainous areas, few field observations exist to supply ground referenced values. This lack of ground observations is especially true at higher latitudes (Fagre et al., 2003; Tucker et al., 2001). However, previous studies have found that topographic variables play a key role in modeling environmental processes such as snowmelt and accumulation; vegetation spatial patterns and life cycles; and available energy through solar radiation (Bacher & Jeanneret, 1994; Hopkins, 1918; Walker et al., 1995). Here we analyze and compare the relationships between start of season and physical, topographic, and land cover variables including temperature, elevation, aspect, slope and ecosystems including tundra, evergreen, deciduous, and mixed forests. We compared the results with Hopkins' bioclimatic law – which was originally developed for all areas east of the Rockies – and states that start of season (SOS) is delayed 3.3 days for every 100 m increase in elevation, 4 days for every one degree change in latitude northward and 1.25 days for every one degree change in longitude westward (Hopkins, 1918). In addition, for 2007, we examine the effect of monthly temperature on the date of spring green-up. 2. Study area The study area consists of six zones dispersed throughout the central and northern Rocky Mountains of both the United States and Canada (Figs. 1, 7) spanning a latitudinal range between 30°18′27 N and 54°18′42 N. Three of the six zones encompass the interior and boundaries of national parks while the remaining three are located within national forests or ecological research stations. A table is provided with more detailed characteristics (Table 1). 2.1. Jasper National Park (JNP) JNP is the largest of the parks in the Canadian Rockies and spans over 10,878 km2 ranging from 53°28′29 N to 52°8′10 N. The three dominant vegetation zones are dryland forest, subalpine evergreen forest, and alpine tundra. The average elevation throughout the park is 2332 m with the highest mountain peak at 3684 m. The average yearly temperature is 3.4 °C; total annual precipitation is 31.9 cm (National Climate Archive, 2007). 2.2. Glacier National Park (GNP) GNP is located at approximately 48°N with an elevation gradient extending from 800 m to over 3000 m. The park spans an area of 4058 km2 and is split by the Continental Divide. The west side of the park is influenced by cool moist Pacific maritime weather whereas east of the Continental Divide weather conditions are driven more by dry windy continental air masses arriving from the Canadian provinces. Dominant vegetation zones are evergreen forest on the mountain slopes, and alpine tundra. GNP receives an average annual precipitation of 121 cm and has an average annual temperature of 3.9 °C (PRISM, 2010). 2.3. Lower Lewis Range (LLR) The LLR begins just south of GNP and covers approximately 12,226 km2 between 48°15′55 N and 47°8′58 N. Like GNP this area is split by the Continental Divide and experiences an average elevation of 814 m with the highest point at 1685 m. An average temperature of 4.9 °C and total annual precipitation of 88 cm support grasslands, shrublands, evergreen forest, and alpine tundra (PRISM, 2010). 2.4. Bitterroot Range (BR) This range is located in the panhandle of Idaho and westernmost Montana and is dominated by grassland and shrublands, montane

1222

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

Fig. 1. Six study zones selected in the western United States and southwestern Canada. These sites span a range of latitudes from 53°28′29 N to 39°44′8 N; elevations of 897 m to 3971 m, and slopes as little as 0° to as steep as 88.32°. SNOTEL sites analyzed are indicated by the gray triangles. The full mosaiked MODIS boundary is designated by the light gray dashed line.

forest, and alpine tundra. The study zone encompasses an area of 3486 km² beginning at 46°38′55 N and running south to 45°51′39 N. The average annual temperature is 5.6 °C and the total annual precipitation is 103 cm (PRISM, 2010).

higher than other sites at 2524 m with a maximum elevation of 4195 m. Dominant vegetation zones throughout the park are grass and sagebrush valleys, montane forest, and alpine tundra.

2.5. Yellowstone National Park (YNP)

2.6. Niwot Ridge Long Term Ecological Research Station (NR)

YNP is the core of the Greater Yellowstone Ecosystem and covers 8987 km², mostly in the northwestern corner of Wyoming beginning at 45°N.The average annual precipitation throughout YNP is 77 cm while the average temperature is 2.9 °C (PRISM, 2010). Elevation changes are less drastic in this zone; however the average elevation is

The Niwot Ridge Long Term Ecological Research Station covers an area of 3397 km2; has an average elevation of 2903 m and is home to extensive alpine tundra, montane forest, subalpine forest, and open meadow. NR and the surrounding area within the study zone (40°22′ 41 N–39°44′8 N) are characterized by an average temperature of 3.9 °C

Table 1 Common characteristics of the six chosen study zones. Location

Size (km2)

Location

Avg. elev (m)

Avg. temp (°C)

1. Jasper NP

10,878

53°28′29 N–52°8′10 N

2332

3.4

2. 3. 4. 5. 6.

4058 12,226 3486 8987 3397

48°59′51 N–48°14′51 N 48°15′55 N–47°8′58 N 46°38′55 N–45°51′39 N 45°1′49 N–44°8′17 N 40°22′41 N–39°44′8 N

1568 814 1821 2524 2903

3.9 4.9 5.6 2.9 3.9

Glacier NP Lower Lewis Range Bitterroot Range Yellowstone NP Niwot Ridge LTER

Total precip (cm) 31.9 121 88 103 77 64

Dominant landcover Alpine tundra, dryland forest, and subalpine evergreen forest Alpine tundra and evergreen forest Alpine tundra, evergreen forest, grassland, and shrubland Alpine tundra, montane forest and shrubland Alpine tundra, montane forest, grass and sage brush valleys Alpine tundra, subalpine forest and montane forest

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

throughout the year, increased solar radiation, and lower annual precipitation than the other study zones at 64 cm (PRISM, 2010). 3. Datasets and processing 3.1. MODIS data The MODIS instrument provides near-daily repeated coverage of the earth's surface with 36 spectral bands and a swath width of approximately 2330 km. Seven bands are specifically designed for land remote sensing with a spatial resolution of 250 m (bands 1–2) and 500 m (bands 3–7). Each MODIS swath is broken into 10 by 10 degree tiles that are numbered vertically and horizontally. For this study we selected the Nadir BRDF-Adjusted Reflectance (NBAR) data set with a spatial resolution of 500 m (MCD43A4v5) covering the ISIN tiles h9v04, h9v05, h10v03, and h10v04. The MCD43A4v5 product is created with the use of a bidirectional reflectance distribution function which models values so as to represent the scene from a nadir view (Lucht, 1998; Lucht et al., 2000; Lucht & Lewis, 2000; Lucht & Roujean, 2000; Lucht et al., 2000). MODIS NBAR is based on multidate, cloud-cleared and atmospherically corrected surface reflectance (MOD09) acquired over a 16-day period. NBAR values are computed for each of the seven spectral bands and normalized to the mean solar zenith angle of each 16-day period. Since view angle effects and cloud and aerosol contamination have been minimized, MCD43A4 delivers a more stable and consistent dataset better suited for global land surface modeling (Lucht et al., 2000; Lucht et al., 2000). This is of particular value for high latitudes and mountainous areas where the availability of high quality imagery is constrained by cloud cover and extreme variations in sun-target-sensor geometry. Data was downloaded for the dates between January 1st 2007 and December 31st 2007 including 46 composites at an eight day temporal resolution. We used the binary MODIS quality flags to select only cloud-free data of ideal quality, and calculated the normalized difference infrared index (NDII) using the NIR (band 2 centered at 858 nm) and MIR (band 7 centered at 2130 nm) (de Beurs & Townsend, 2008; Hardisky et al., 1983; Hunt & Rock, 1989) and normalized difference vegetation index (NDVI) calculated using the NIR and the red reflectances (band 1 centered at 645 nm). The calculations for each were performed as follows: NDII =

Band2−Band7 Band2 + Band7

NDVI =

Band2−Band1 : Band2 + Band1

ð1Þ

Both of these indices allow for the observation of vegetation dynamics and health based on levels of chlorophyll or water content in the foliage. After calculating the indices, we resampled the data to 16-day temporal resolution by selecting maximum NDVI/NDII values. Thus, a total of 23 scenes were used for the land surface phenology analysis. We re-projected the data into geographic coordinates and mosaicked based on the georeferenced location to create one continuous scene between latitude 30°18′27 N and 54°18′42 N and longitude 99°49′31 W and 120°2′40 W (Figs. 1, 7).

1223

International Geosphere Biosphere Programme (IGBP) classification scheme, which includes 11 natural vegetation classes, three developed and mosaiked land classes, and three non-vegetated land classes. Band10 represents the secondary classification schemes based on the IGBP data. Land cover bands were stacked with the DEM derived outputs for elevation, slope, aspect, and solar radiation as well as the SOS season layer, and four vegetation classes were chosen based on their relevancy and abundance across our study region. These classes include tundra (23%), evergreen forest (43%), deciduous forest (6.3%), and mixed forest (18%). 3.3. SNOTEL data In situ measurements of daily historical average temperature data were obtained from the automated SNOTEL (snow telemetry) weather station network (United States Department of Agriculture, National Resources Conservation Service (http://www.wcc.nrcs.usda. gov/snotel/)). There are over 600 SNOTEL stations in 13 states primarily in the western area of the United States. SNOTEL ground based weather stations measure climatic data including daily minimum and maximum temperature, as well as snowpack parameters such as snow water equivalent, snow density, and snow depth. As SNOTEL sites currently only exist within the U.S., we collected daily minimum and maximum temperature for five of the six study areas used in this study. Within these states, 80 sites were sampled for Montana, 81 sites in Idaho, 80 sites in Wyoming, and 86 sites in Colorado (Fig. 1) for the months of March through July. As site locations were limited within the individual study zones (Fig. 1), we chose to examine all stations by state. Average daily temperature based on the minimum and maximum temperatures was calculated daily for each station. Monthly averages were derived by taking the mean of this result. 3.4. Digital elevation model In order to examine the relationship between the topographic variables and the timing of green-up or start of season (SOS) in mountain vegetation, two elevation datasets were acquired. The first elevation model was obtained from the United States Geological Survey's (USGS) Earth Resources Observation and Science (EROS) data center's National Map Seamless server (http://seamless.usgs. gov/index.php). This USGS National Elevation Dataset (NED) was retrieved at a spatial resolution of 30 m for Idaho, Montana, Wyoming, and Colorado. As the coverage of the NED does not extend into Canada, a second 90 m dataset, The Shuttle Radar Topography Mission (SRTM) DEM, was obtained from the USGS Earth Explorer server (http://edcsns17.cr.usgs.gov/EarthExplorer/) for the area encompassing Jasper National Park. These two elevation products were aggregated and scaled up to allow for the calculation of slope, aspect, and solar radiation at a spatial resolution corresponding to the scale of the NDVI data. The final outputs were clipped to cover the appropriate study area. 4. Methodologies

3.2. MODIS land cover data 4.1. Slope, solar radiation, and aspect To examine the possible influence of varying land cover types on SOS we used the MCD12Q1 land cover product. This product was obtained at a spatial resolution of 500 m for the entire western mountain range. The MCD12Q1 product contains 16 bands each representing a different classification scheme describing land cover properties of annual observations obtained from both the Terra and Aqua satellites (Friedl & Brodley, 1997; Friedl et al., 2002; Muchoney et al., 2000). All tiles were re-projected to geographic coordinates and mosaicked together creating one seamless image. Bands 1 and 10 specifically were used for the analysis. The first band represents the

Both digital elevation models consisted of a regular matrix of elevation values from which topographic functions including slope, aspect and solar insolation were calculated within a GIS and rendered for visualization and analysis. Each function's output was produced as an individual layer of values and was joined with the output for SOS as discussed below. The slope and aspect were calculated using the ArcGIS Surface tool within the Spatial Analyst Toolset (ESRI, 2006). Within the GIS, slope identifies the rate of maximum change in altitude, from each cell. The

1224

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

mask the increase due to green-up. This will lead to a derived SOS value later than the actual SOS. Tgreening = max ðt€½0; 200ÞðNDIIðtÞbNDIImin + εÞÞ:

ð3Þ

Here ε is chosen as follows: 1) it must be larger than the noise affecting the NDII time profile, and 2) it must be smaller than the first NDII increase due to vegetation growth. In addition, we insert a value of 200 as we assume that SOS will be prior to this day (July 18th) (Delbart et al., 2005). 4.3. Regression analysis

Fig. 2. NDVI and NDII for 2007. This pixel was located in Yellowstone National Park 44°54′30 N and 111°2′40 W. NDII values are highest when snow is on the ground and begin to decrease as snowmelt occurs. In contrast, NDVI values are lowest when snow cover is present and begin to increase when snowmelt occurs and vegetation green-up begins. SOS for NDII is estimated as the date that the NDII signal begins to increase again. NDVI SOS is estimated as the date when the NDVI signal has reached fifty percent of its maximum value.

output slope values are in percent slope to allow for their inclusion in a linear function. The values of 0% slope therefore correspond to flat areas while values of 100% correspond to an area that is completely vertical. Aspect identifies the downslope direction of the maximum change in altitude from each pixel to its neighbors. The output values were expressed in positive azimuth degrees from 0° to 359.9°, and were measured clockwise from the north. For our analysis aspect was transformed to a linear variable using the following cosine function (Beers et al., 1956): Aspectlin = cosð45−AÞ + 1

ð2Þ

where aspectlin equal to −1 indicates northern slope, and aspectlin equal to 1 represents a southern slope. To derive insolation we used the solar radiation tool within the ArcGIS surface analysis toolset (Rich et al., 1994). Solar radiation was calculated based on the area of insolation across each study zone. Calculations were performed for the entire year of 2007 at monthly intervals. 4.2. Land surface phenology analysis Over the last few years, several methods have been developed to detect the onset of green-up from vegetation index data (Reed et al., 1994; White et al., 1997, 2009; Zhang et al., 2003). In order to characterize the land surface phenology here we calculate start of season based on both NDVI and NDII (Fig. 2). Comparative studies have shown that thresholds based on variations of an NDVIratio consistently estimate start of season (SOS, in days) that corroborates ground-based phenology data (White et al., 1997, 2009). Thus, for the determination of SOS based on NDVI we used the Midpointpixel method, where the value midpoint was set at fifty percent of the maximum NDVI (White et al., 1997). SOS is assumed to have occurred on the date that the ratio of NDVI rises past the 50% point. Since NDVI increases with snow melt, SOS in snowy areas based on the Midpointpixel method is occasionally estimated earlier than actual vegetation green-up. With the addition of the middle infrared (MIR) band of the MODIS sensor, potential snow effects on the estimation of SOS should be better controlled. Snow melt can be seen as a decrease in NDII values and subsequent vegetation green-up results in an increase in values (Fig. 2) (Delbart et al., 2005). In situations where vegetation green-up occurs during snowmelt the NDII decrease can

To relate SOS to the physical parameters present within our study zones the SOS output layer was overlaid with the DEM and the GIS derived outputs for slope and solar radiation. A random sample of ten percent of the total number of pixels was selected for each study zone. In addition to the above variables, longitude and latitude were extracted for each randomly selected pixel. Since latitude and longitude changes are relatively minor in our study areas (less than 15° from North to South and East to West), we assume that a linearization of the geographic coordinates is not necessary before including them into our linear regression models. A backward stepwise linear regression model (Eq. 4) was used to determine the independent variables with the highest coefficients of determination: SOS = intercept + αelevation + βlatitude + γlongitude + δslope + ζsolar radiation:

ð4Þ

In addition, we applied a backward stepwise linear regression to examine predicted and derived SOS for four aspects independently: north, south, east and west. North facing slopes were designated as any aspect value falling between 315° and 45°, easterly slopes were between 46° and 135°, southerly slopes were between 136° and 225°, and westerly slopes were between 226° and 314°. This analysis was performed for the entire western mountain range. The next step examined the relationship between temperature and SOS. For this step SNOTEL temperature data were used. The model included monthly average temperatures for March, April, May, June, and July, as well as elevation, slope, longitude and latitude. In the final step, we investigated the correlations between vegetation and SOS. This dataset was sampled representing the entire western mountain range using the exact samples and lines as in the NDVI and NDII non stratified analyses. If the dominant class of a pixel fell outside our classes of interest (Tundra, Deciduous Forests, Evergreen Forests, and Mixed Forests), we investigated the second-dominant class depicted in band 10. The addition of the secondary land cover classes allowed for an increased diversity in the samples and especially boosted the percentage of deciduous and mixed forests. The final model included elevation, slope, aspect, longitude, latitude, and the relevant land cover. 5. Results 5.1. Individual study zone SOS Table 2 shows the relationship between the satellite derived SOS from NDVI (Fig. 3) and elevation(m), slope (%), aspectlin, and spatial location expressed as longitude and latitude for our six study zones. Results indicate that SOS can be determined for all zones independently with a significant coefficient of determination (R2adj values between 0.32 and 0.69). The best fitting model, with a R2adj of 0.69, was found in NR, the study area farthest south. NR is also the smallest study area with a sample size of 956 points. Despite the southerly location of this site, the average elevation is consistently higher than at most other more northerly locations and thus there are SOS dates

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

1225

Fig. 3. Predicted SOS date based on significant variables, against actual SOS as determined using the NDVIratio method. Study zones run north to south. The vertical black line represents the 1:1 relationship between the predicted and actual SOS values. The SOS can be predicted for all the zones with a significant adjusted R2.

later than in LLR, BR, and YNP. JNP, the most northern of the six study areas, has the latest SOS dates. GNP had the lowest R2adj at 0.32. The relationship between elevation and SOS within all six study areas was consistent with the relationship between elevation and start of season as recorded in the eastern deciduous forests by Hopkins' Bioclimatic Law (1918). This law states that start of season is delayed, 3.3 days for every 100 m increase in elevation. We found a SOS delay for every 100 m increase in elevation ranging between 2.9(±0.5) days in GNP and 6.8(±0.1) days in NR. Relationships found between SOS and latitude were significant for JNP, LLR, BR, YNP, and NR, with LLR and BR indicating a delay in SOS between 12(±1.9) and 14(±0.8) days with a change of one degree southward and a delay of 23(±0.7) days for every one degree shift northward in YNP. The findings for the relationship between SOS and longitude show that all parks north of Wyoming experience a delay in SOS between 5.8(±1.3) and 53 (±3.4) days with a change of one degree westward. SOS in JNP, GNP,

LLR, and YNP exhibited significant correlations with aspect. In JNP, LLR, and YNP a movement of 1 rad or 57° from the southeastern aspect, results in a delay of SOS between 1(±0.2) and 5(±0.3) days. In contrast we find that in GNP the same shift in aspect results in a 1.8 (±0.4) advance in SOS. Table 3 and Fig. 4 show the relationship between satellite derived SOS from NDII and the model estimates based on longitude, latitude, elevation, slope and aspect. The relationship between elevation and SOS also follows closely with Hopkins' Bioclimatic Law showing a delay ranging between 3.2(±0.7) and 6(±0.1) days for every 100 m increase in elevation. As with NDVI derived SOS, the pattern shown for both longitude and latitude does not follow closely with the findings of Hopkins, with the exception of GNP and LLR where a shift in one degree northward causes a delay in SOS between 8(±1.7) and 10 (±0.7) days. In JNP and YNP the relationship is reversed so that a one degree change in latitude southward results in a SOS delay on average

1226

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

Fig. 4. Predicted SOS date for all six individual study zones based on significant variables, against actual SOS as determined using NDII data. The vertical black line represents the 1:1 relationship between the predicted and actual SOS values.

of 15 days. A shift in longitude of one degree westward will delay SOS in JNP, LLR, BR, YNP, and NR between 11(±0.7) and 30(±3) days. We also found that an eastern or western change of 57° in the southeastern aspect resulted in a delay in SOS between 1(±0.7) and 4(±0.2) days. 5.2. Western Mountain Range SOS Results (Tables 2, 3) show that NDVI derived SOS dates can be predicted for the entire western mountain range with a significant R2adj of 0.56 (Fig. 5) and for NDII derived SOS with a slightly higher R²adj of 0.58 (Fig. 6). We found that NDVI SOS and NDII SOS do not exhibit significant differences in the predicted delay of spring greenup based on an equal positive change in elevation. NDVI SOS is delayed 3.6(±0.2) days for every 100 m increase in elevation while NDII SOS outputs indicate a delay of 3.9(±0.3) days for the same increase in elevation. Longitude was also significant for both NDVI and

NDII indicating a SOS delay of approximately seven days (NDVI), and 6 days (NDII) for every degree westward, latitude was significant for both NDVI and NDII SOS. For NDVI we found SOS delay of 1.7(±0.2) days for every degree in latitude southward, while for the NDII SOS, we found a delay in SOS of 1.3(± 0.1) days for every degree northward. The results for NDII SOS follow with Hopkins' bioclimatic law. Finally, we find that for both NDII and NDVI derived SOS a change to the east or west in southeastern aspect of 57° causes a 2.4(±0.4) day delay in SOS. To further understand the effect of aspect, we divided the NDVI SOS data into four categories, North (315°–45°), East (45°–135°), South (135°–225°) and West (225°–315°). The results (Table 4) from this additional analysis reveal the relationship between NDVI SOS, and parameters including elevation, longitude and latitude vary by aspect. In addition, the strength of the relationship also varies (R²adj values range between 0.35 and 0.50) and is strongest for south facing slopes and weakest for north facing slopes. We found that for every 100 m

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

1227

forest when examining the NDVI based SOS. Aspect was relevant in both analyses across all vegetation types except the mixed forest. For both NDVI and NDII we found that a shift of 1° aspect results in a delay of SOS between 2(± 0.1) and 3(±0.5) days. Slope was significant in both analyses, however for the NDVI analysis it is only relevant for tundra and mixed forest while the NDII data shows relationships among tundra and evergreen (Figs. 9, 10). 5.3. State-wide analysis at SNOTEL locations

Fig. 5. Results derived from a linear regression model for the entire western mountain range based on the NDVI derived SOS. The vertical black line represents the 1:1 relationship between the predicted and actual SOS values.

increase in elevation on a northerly slope, SOS is delayed 5.2(±0.2) days, while such an elevation increase on south facing slopes results in a delay in SOS of 6.1(±0.2) days. Latitude was significant for all zones and a one degree shift northward delays SOS between 4.3(±0.1) days to 6.5(±0.5) days. Longitude was significant for all slopes except the north facing, and a one degree shift westward delays SOS between 5.1 (±0.7) and 7.1(±0.6) days. On average, SOS was as much as 20 days earlier on south facing slopes than those slopes on the north side of the mountains. To study the relationship between varying vegetation zones and SOS our data was stratified based on four land cover classes; tundra, evergreen forest, deciduous forest, and mixed forest. The results for both the NDVI and NDII analyses (Tables 6, 7) indicate that elevation is the most significant driver of SOS, despite the predominant land cover present, and consistently follows Hopkins' law (a shift in 100 m north results in a delay between 3(±0.2) and 4(±0.1) for both the NDVI and NDII analyses). When examining the SOS derived from NDVI data, we found the highest model fit (R²adj = 0.62) in deciduous forests. The lowest model fit was found for evergreen forests (R²adj = 0.37). For the NDII output, the tundra areas revealed the highest model fit (R²adj = 0.63), while the mixed forest revealed the lowest model fits with an R² of 0.52. For both the NDVI and NDII analyses longitude proved to be significant across all land covers with shifts in SOS between 5(± 0.1) and 8(±0.4) days for the NDVI data and 4(±0.2) and 6(±0.6) days for the NDII data. Latitude proved to be the least significant variable only exhibiting correlations among the mixed

Fig. 6. Results derived from a linear regression model for the entire western mountain range based on the NDII derived SOS. The vertical black line represents the 1:1 relationship between the predicted and actual SOS values.

Parameter estimates from the backward stepwise linear regression (Table 5) show that NDVI SOS is strongly affected by a change in elevation for the SNOTEL locations available in Montana, Idaho, Wyoming, and Colorado. Delays in SOS based on elevation are similar to those seen in the NDVI and NDII analyses (Tables 2, 3). The relationship with latitude indicates that SOS was delayed the least in Colorado 3.3(±0.2) for every degree northward and the most, 8.6 (± 0.4) days, in Montana. Monthly average temperature is a significant factor in determining SOS in Montana, Idaho, and Wyoming; however the months that provide the most significant input between all three states vary and reveal an interesting pattern. Montana, with SOS dates ranging from April 1st through June 3rd, shows that an increase in monthly average June temperature of one degree Celsius initiates SOS about 17(±0.4) days earlier, while the same increase in the average temperature in May will delay SOS about 17(±0.3) days. Idaho exhibits SOS dates as early as March 21st and no later than May 25th and the most significant correlation with temperature appears to be in April where an increase of one degree Celsius leads to an earlier SOS of 10.3(±0.1) days. Interestingly, an increase of one degree Celsius in March delays SOS by 8.7(±0.2) days. Wyoming's SOS dates range from April 9th to May 27th. SOS dates throughout Wyoming show similar associations with the average temperatures in April and May. An increase of one degree Celsius in April delays SOS 6.3(±0.6) days while a one degree increase in May results in an advance in the SOS of 6.4(±0.3) days (Fig. 8). 6. Discussion 6.1. The use of NDVI and NDII for determining SOS With the use of data collected via satellite sensors, vegetation indices can be produced that allow for the monitoring and characterizing of vegetation phenology patterns in remote areas such as mountainous environments. However, characterization of land surface phenology in mountainous environments is hindered by varying topography, satellite angle when imagery was collected, and a lack of ground-based observations for validation. We have shown that both NDVI and NDII data can be used to predict SOS with a significant R2 for each study zone despite the variability present in land cover, elevation shifts, and general topographic characteristics. We found that NDVI works best for predicting later SOS values which are likely to occur at higher elevations and NDII is more useful when predicting earlier SOS which we assume to be associated with pixels sampled at lower elevations. In several areas we found that NDVI estimates SOS earlier than we would expect based on the physical variables that we included in our study. We assume that this discrepancy results from the premature increase in NDVI as a result of snowmelt (Delbart et al., 2005, 2006). In the context of boreal regions, Delbart et al. (2005, 2006) found that NDWI, which is comparable to NDII, was more efficient in estimating the date of onset of greening than methods based on NDVI. However, the use of the NDII method can lead to a delay in SOS in areas where snow cover is prevalent later in the year and vegetation green-up may actually coincide with the snowmelt. At lower elevations, snowmelt is more likely to occur earlier which allows stability to remain in the NDII signal.

1228

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

Fig. 7. Depictions of the satellite derived SOS (in Julian days) across our study zone based on (A.) NDVI and (B.) NDII. The areas in white indicate no data. The SOS is said to have occurred later as green darkens across the area.

6.2. The effect of elevation on SOS SOS in vegetation is controlled by several forces including photoperiod, temperature, and precipitation (Bacher & Jeanneret 1994; Hopkins, 1918; Walker et al., 1995). In mountain environments, shifts in elevation and aspect can occur over very short distances altering local temperature patterns. In addition, shifts in latitude and longitude across our study area affect both photoperiod and precipitation. These variations contribute to the specific patterns and processes among the vegetation present not only in mountain environments but in the landscapes bordering these mountains (Bacher & Jeanneret 1994; Barry, 1992). When examining the zones independent of one another and without the inclusion of temperature, we found that elevation was consistently the most significant variable

and was highly correlated to SOS as derived from both the NDVI and the NDII data. This result corresponds to the findings from Hopkins (1918) which stated that elevational changes proved to be the most significant determinant in seasonal variations within vegetation. Similar relations were also found when the results were examined at a regional scale (the entire range of the study zones) and the state-wide scale which included temperature measurements. 6.3. The effect of aspect The aspect of a slope can have very significant impacts on vegetation type, soil type, snowmelt, insolation, runoff, and fire patterns, which may all lead to variations in the timing of phenological events (Bliss, 1956; Mark, 1970; Ratcliffe & Turkington, 1989). In the northern hemisphere, south facing slopes receive more direct sun and are subsequently drier while north facing slopes experience lower levels of evapotranspiration due to decreased sunlight and cooler winds. These cooler winds lead to microclimates that experience increased snow accumulation and later snowmelt (Barry, 1992; Wagner & Reicheggar, 1997; Walsh, 2004). Our results found that for the south facing slopes the average SOS was as much as 3 weeks earlier than the SOS on northern slopes (Table 4). This is not an unexpected finding as the consistently higher temperatures on south facing slopes as well as the lower levels of snow accumulation would encourage earlier vegetation growth. 6.4. The effect of temperature on SOS

Fig. 8. Results derived from a linear regression model for all four states at the SNOTEL locations. Observed SOS values are based on the NDVI SOS. The vertical black line represents the 1:1 relationship between the predicted and actual SOS values.

The relationship between temperature and vegetation phenology is of increased interest as researchers continue to analyze the biotic implications of responses to climate warming (Schwartz et al., 2006; Sphen & Korner, 2002). These responses are of particular interest in high elevation and high latitude environments due to recorded recent warming trends and anticipated continued warming, (Myneni et al., 1997; Schwartz et al., 2006). Most organisms are adapted to live within a relatively narrow range of temperatures and will not thrive if temperatures breach this threshold. The boundaries imposed by

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

1229

Fig. 9. Predicted SOS date based on significant variables, against actual SOS as determined using the NDVIratio method. Samples are stratified based on the four dominant land covers across the entire western mountain range including Tundra, Evergreen, Deciduous, and Mixed Forest. The vertical black line represents the 1:1 relationship between the predicted and actual SOS values.

temperature on the growing season of plants, is an important example of this dependent relationship. For this analysis we found that when temperature was examined along with the additional proxy variables, elevation remained the most significant predictor raising the R2adj values as much as 0.31 points when included in the regression. This is most likely a result of different vegetation types at different elevations. However, even with the inclusion of elevation, significant relationships were found between SOS and monthly average temperatures (R²adj values of 0.45 were found when examining SOS against temperature and R²adj values as high as 0.68 when examining the relationship between elevation and temperature). Existing research has shown that higher temperatures occurring either at the onset of green-up or early in the growing season can speed plant development and contribute not only to a more rapid onset of spring, but also to a lengthened growing season (Menzel & Fabian, 1999; Myneni et al., 1997; Peñuelas & Filella, 2001). For Montana, Idaho, and Wyoming results indicated that an increase of one degree Celsius, during the early growing season, advances SOS anywhere from 16.9 days (Montana) to 6.4 days (Wyoming). In contrast, when temperatures increase in the months that precede SOS, spring can be delayed up to nine days. We suggest that this pattern may be attributed to one of two very different weather events — a late frost, or increased winter snow fall. An underlying hypothesis found in the literature (Augspurger, 2009; Cannel & Smith, 1986; Gu et al., 2007, 2008; Meehl et al., 2000) reports on the deleterious effect of late frost events on premature vegetation growth. These frost events occur when mild winter temperatures lead to warm early springs, brought on by a warming climate. The warmer temperatures initiate an earlier green-up and the late frost leads to a subsequent loss of venerable new vegetation. This initial loss to vegetation development can delay the recorded SOS in both ground and satellite observations as plants strive to regain stamina in their physiological processes (Gu et al.,

2008). The second hypothesis relates to a type of Clausius–Clapeyron relationship which shows that higher atmospheric temperatures lead to an increase in atmospheric moisture. When the temperatures cool, this moisture is released as snow, thus resulting in a positive relationship between air temperature and snow (Davis et al., 1999). An additional study by Changnon et al. (2006) reported that based on snowfall and temperature measurements between the years of 1901 and 2000, the most extreme snowstorms have occurred in warmer than normal years. With increased or heavier snowfall occurring in months preceding SOS, vegetation green-up is likely to be delayed. Another explanation for the earlier and later SOS patterns within our results may be the varying vegetation classes among samples. In Yellowstone National Park alone there are six dominant vegetation classes including both coniferous and deciduous forests as well as grass and sage brush valleys. Seasonal green-up patterns will differ dependent on these vegetation classes, some experiencing much earlier SOS (grass and sage brush), while others experience much later SOS (alpine tundra and montane coniferous forest) (Bradley & Mustard, 2008). The lack of correlation found in Colorado between SOS and temperature may be attributed to the uniformly high elevational range of SNOTEL locations. In Colorado, the average elevation is 3126 m which is significantly higher than that found in Montana (2126 m), Idaho (1992 m), and Wyoming (2643 m). With minor elevational variation, temperature remains relatively constant resulting in negligible impact on SOS (Figs. 8, 9). 6.5. SOS for different land cover classes It has been shown that studying trends in plant growth throughout boreal regions is critical in understanding ecosystem functions and feedbacks between land cover and photosynthetic activity. Goetz et al. (2005) and Goetz and Prince (1996) indicated that gross

1230

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

Fig. 10. Predicted SOS date based on significant variables, against actual SOS as determined using NDII data. Samples are stratified based on the four dominant land covers across the entire western mountain range including Tundra, Evergreen, Deciduous, and Mixed Forest. The vertical black line represents the 1:1 relationship between the predicted and actual SOS values.

photosynthetic activity has been increasing since 1982 in high latitude vegetation and that this trend can be attributed to an earlier SOS and subsequent lengthening of the growing season. While this work only examines the initiation of the growing season for the year 2007, it can be seen that the effect of biophysical parameters on the day of spring green-up differs significantly by land cover class. SOS could be estimated with highest significance in the tundra and deciduous forests in both the NDVI and NDII analyses. This finding is not surprising as we would expect both indices to record the strongest increases in spring vegetation growth within these two zones. Evergreen forests are the dominant land cover in these high latitude

Table 2 Parameter estimates for each significant independent variable as it relates to NDVI SOS within the specified study zone. The sample size is an effect of ten percent of the size of the study zone. Negative longitudinal values indicate a delay in SOS westward while negative latitudinal values indicate a delay in SOS southward. In addition to these variables, the amount of solar radiation present did significantly influence SOS predictions for GNP during the months of March, May, and July; LLR for the month of February only; and BR for both February and March. SOS in YNP and NR showed no significance to solar radiation or slope. Missing estimates are not significant. NDVI Elevation (m) Latitude (°) Longitude (°) Slope (%) Aspect N Adj. R² p-value

JNP

GNP

LLR

BR

YNP

NR

All zones

0.032 − 23.75 − 21.17 N/A − 3.62 9253 0.51 b 0.001

0.029 N/A − 5.8 N/A 3.13 2611 0.32 b0.001

0.043 14.48 − 31.40 − 0.15 − 1.77 5184 0.41 b 0.001

0.060 12.12 − 52.98 − 0.09 N/A 1113 0.49 b 0.001

0.052 − 23.19 N/A − 0.39 N/A 4108 0.58 b0.001

0.068 7.34 N/A N/A 1.95 956 0.69 b 0.001

0.036 − 1.73 − 6.90 N/A − 2.42 23,225 0.56 b 0.001

regions comprising 43% of our samples. It is most likely that when examining these pixels we are seeing the initial understory growth appearing as a rapid increase in the signal. As the over story begins to fill in, the increase in the spectral signature slows. It can then be assumed that any remaining incline is attributed directly to the increased density of the arboreal environment (Figs. 9, 10). 7. Conclusion Mountain environments present unique challenges to the characterization of vegetation phenology because of the extensive variability in elevation, temperature, slope, aspect, solar radiation, and land cover. Observations derived from data collected onboard space born satellite sensors provide excellent opportunities to characterize the land surface

Table 3 Parameter estimates for each significant independent variable as it relates to SOS as determined from NDII data within the chosen study zone. The SOS is delayed 21 days for every degree westward for LLR and 15 days for YNP. In addition to the above parameters, SOS is significantly influenced by slope and solar radiation for the months of February through July only within the LLR. NDII

JNP

GNP

LLR

BR

YNP

NR

Elevation (m) 0.032 0.043 0.056 0.060 0.053 0.049 Latitude (°) − 16.62 10.32 7.94 NA − 15.39 10.67 Longitude (°) − 10.84 NA − 19.70 − 23.93 − 15.03 30.39 Slope (%) NA 0.11 NA NA NA NA Aspect − 3.64 NA − 2.69 NA − 1.44 NA N 9253 2611 5184 1113 4108 956 Adj. R² 0.40 0.43 0.54 0.57 0.38 0.57 p-value b 0.001 b 0.001 b 0.001 b 0.001 b 0.001 b 0.001

All zones 0.039 − 1.25 − 6.22 NA − 2.45 23,225 0.58 b 0.001

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233 Table 4 The average SOS for North, East, South, and West aspects along with parameter estimates for each significant independent variable as it relates to SOS and aspect direction. As expected the SOS values were earlier on south and west slopes on average. The negative longitudinal values indicate a delay in SOS westward while the negative latitudinal values indicate a delay in SOS southward. The amount of solar radiation present did significantly influence SOS predictions for south facing slopes during the months of February, May, and June and for west facing slopes during the months of February and March. Missing estimates are not significant. Aspect

North

East

South

West

Avg. SOS (DOY) Elevation (m) Latitude (°) Longitude (°) N Adj. R² p-value

135 0.052 6.51 N/A 732 0.35 b 0.001

137 0.060 4.46 − 5.14 1744 0.45 b0.001

114 0.061 4.32 − 7.09 685 0.50 b0.001

120 0.055 4.82 − 5.71 1718 0.47 b 0.001

Table 5 Parameter estimates for the SNOTEL temperature analysis indicate that elevation, monthly average temperature, and latitude are strongly significant determinants of SOS at SNOTEL locations in all four states. Based on estimates for elevation, SOS is delayed as little as 3 days in Colorado for every 100 m increase in elevation to as much as 60 days in Montana. Neither slope nor solar radiation significantly influences SOS. SNOTEL

Montana

Idaho

Wyoming

Colorado

All zones

March April May June Elevation (m) Latitude (°) N Adj. R² p-value

N/A N/A 16.94 − 16.95 0.062 8.55 80 0.47 b0.001

8.66 − 10.30 N/A N/A 0.035 4.17 81 0.43 b 0.001

N/A 6.31 − 6.40 N/A 0.041 4.95 80 0.47 b 0.001

N/A N/A N/A N/A 0.048 3.28 86 0.30 b 0.001

1.39 N/A − 2.41 N/A 0.036 3.18 317 0.44 b 0.001

phenology of remote environments worldwide (Delbart et al., 2005; White et al., 1997, 2009; Zhang et al., 2004). Here we compared the SOS estimates based on NDVI, the most commonly used index in studies examining vegetation dynamics, and NDII, another spectral index which has recently been used to reduce the effects of snowmelt on the estimation of SOS. The results of our analysis indicate that both NDVI and NDII can be used with varying degrees of confidence to predict SOS in a North American western mountain range. However, in future studies we aim to develop a methodology which combines properties of both indices. This combination will allow for a more accurate estimation of SOS values in boreal regions which experience large amounts of snowfall and subsequent snowmelt. Future analysis of mountain phenology is needed to enhance our current understanding of the effects that elevation, slope, aspect, solar radiation, temperature, and land cover have on the SOS. Here we found that the seasonal green-up pattern of the western mountain range closely follows the relationship stated by Hopkins Bioclimatic

Table 6 Parameter estimates for each significant independent variable as it relates to NDVI SOS across varying land covers. Sample size is an effect of ten percent of the size of the study zone. Missing estimates are not significant. NDVI

Tundra

Evergreen

Deciduous

Mixed

Elevation (m) Latitude (°) Longitude (°) Slope (%) Aspect N Adj. R² p-value

0.040 N/A − 7.52 − 0.15 − 2.73 4852 0.59 b0.001

0.034 N/A − 6.47 N/A − 2.41 11,162 0.37 b0.001

0.036 N/A − 6.49 N/A − 2.59 1338 0.62 b0.001

0.033 1.19 − 4.83 − 0.17 N/A 3723 0.55 b 0.001

1231

Table 7 Parameter estimates for each significant independent variable as it relates to NDII across varying land cover. Sample size is an effect of ten percent of the size of the study zone. Missing estimates are not significant. NDII

Tundra

Evergreen

Deciduous

Mixed

Elevation (m) Latitude (°) Longitude (°) Slope (%) Aspect N Adj. R² p-value

0.044 N/A − 5.90 − 0.19 − 2.97 4852 0.63 b0.001

0.040 N/A -4.15 0.08 -2.21 11,162 0.56 b0.001

0.039 N/A − 4.21 N/A − 3.45 1338 0.59 b 0.001

0.040 N/A − 4.32 N/A N/A 3723 0.52 b 0.001

Law (Hopkins, 1918) in regard to elevation and start of season (at both the local and regional scales and whether accounting for aspect or not) as well as latitude and SOS (at the regional scale only). Aspect proved to be relevant for most parks at the local scale with both NDVI and NDII SOS, as well as for the entire western mountain range. While we did not find consistent significant relationships between SOS and slope, and no direct relationship between SOS and solar radiation, it has been shown that these factors can strongly affect vegetation characteristics including phenology (Malingreau, 1986; Moulin et al., 1997; Schwartz & Reiter, 2000). When studying the SOS within the four dominant land cover classes we found that all four exhibit significant relationships between the measured biophysical variables and the timing of vegetation green-up, however the level of correlation varies depending on the cover type. SOS can be predicted with the highest R²adj for Tundra and Deciduous forests. With the use of finer resolution satellite data it may be possible to ascertain patterns and trends not revealed at a resolution of 500 m. The new Web-Enabled Landsat Data (WELD) project fuses the Landsat Thematic Mapper data with MODIS land products to systematically generate a seamless mosaic dataset at monthly, seasonal, and annual temporal resolutions (Roy et al., 2009). The resulting output is a high spatial resolution product (30 m), with 7 bands which can be used to create numerous spectral indices for vegetation phenology studies. The availability of such a dataset provides an opportunity to study patterns and processes at a higher spatial resolution which may allow a representation of the relationships between variables that are influenced by minor changes in spatial location but strongly contribute to processes in vegetation dynamics. Characterizing the land surface phenology of mountain environments is a complex task. This work is a beginning step to the gathering of a complete understanding of the phenological patterns and process occurring in these highly sensitive ecosystems.

Acknowledgements This work was funded by the following sources: The Geography Department of Virginia Tech; The Sidman P. Poole Endowment; and The Virginia Tech Graduate Student Assembly Graduate Research and Development Program Grant. Additional acknowledgements and thanks are extended to P. de Beurs for the application development that allowed us to estimate the land surface phenology parameters more efficiently. Dr. Jim Campbell and Dr. Lynn Resler are also gratefully acknowledged for their advice, edits, and subject matter contribution.

References Allen, T. R., & Walsh, S. J. (1996). Spatial and compositional structure of the alpine treeline ecotones. Photogrammetric Engineering and Remote Sensing, 62, 1261−1268. Augspurger, C. K. (2009). Spring warmth and frost: phenology, damage and refoliation in a temperate deciduous forest. Functional Ecology, 23, 1031−1039.

1232

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233

Bacher, F., & Jeanneret, F. (1994). Phenology as a tool in topoclimatology. In M. Beniston (Ed.), Mountain environments in changing climates (pp. 217−226). New York: Rutledge Press. Barry, R. (1992). Mountain weather and climate (pp. 402). Cambridge, UK: Cambridge University Press. Beers, T., Dress, P., & Wensel, L. (1956). Aspect transformation in site productivity research. Journal of Forestry, 64, 691−692. Bernbaum, E. (1997). The spiritual and cultural significance of mountains. In B. Messerlie, & J. D. Ives (Eds.), Mountains of the world: A global priority (pp. 39−60). Oxford: Parthenon Publication Group. Bliss, L. C. (1956). A comparison of plant development in microenvironments of arctic and alpine tundras. Ecological Monographs, 26, 303−337. Blyth, S., Groombridge, B., Lysenko, I., Miles, L., & Newton, A. (2002). Mountain watch: Environmental change and sustainable development in mountains. UNEP-World Conservation Monitoring Center. Bradley, B., & Mustard, J. (2008). Comparison of phenology trends by land cover class: A case study in the Great Basin, USA. Global Change Biology, 14, 334−346. Cannel, M. G. R., & Smith, R. I. (1986). Climatic warming. Spring budburst and frost damage on trees. Journal of Applied Ecology, 23, 177−191. Changnon, S. A., Changnon, D., & Karl, T. R. (2006). Temporal and spatial characteristics of snowstorms in the contiguous United States. Journal of Applied Meteorology and Climatology, 45, 1141−1155. Cleland, E. E., Chuine, I., Menzel, A., Mooney, H. A., & Schwartz, M. D. (2007). Shifting plant phenology in response to global change. Trends in Ecology & Evolution, 22, 357−365. Davis, R. E., Lowit, M. B., Knappenberger, P. C., & Legates, D. R. (1999). A climatology of snowfall-temperature relationships in Canada. Journal of Geophysical Research, 104, 11,985−11,994. de Beurs, K. M., & Henebry, G. M. (2004). Land surface phenology, climatic variation, and institutional change: Analyzing agricultural land cover change in Kazakhstan. Remote Sensing of Environment, 89, 497−509. de Beurs, K. M., & Henebry, G. M. (2010a). Spatio-temporal statistical methods for modeling land surface phenology. In I. L. Hudson, & M. R. Keatley (Eds.), Phenological research: methods for environmental and climate change analysis (pp. 177−208). New York: Springer Books. de Beurs, K. M., & Henebry, G. M. (2010b). A land surface phenology assessment of the northern polar regions using MODIS reflectance time series. Canadian Journal of Remote Sensing, 36, S87−S110. de Beurs, K. M., & Townsend, P. A. (2008). Estimating the effect of gypsy moth defoliation using MODIS. Remote Sensing of Environment, 112, 3983−3990. Delbart, N., Kergoat, L., Toan, T. L., Lhermitte, J., & Picard, G. (2005). Determination of phenological dates in boreal regions using normalized difference water index. Remote Sensing of Environment, 97, 26−38. Delbart, N., Toan, T. L., Kergoat, L., & Fedotova, V. (2006). Remote sensing of spring phenology in boreal regions: A free of snow-effect method using NOAA-AVHRR and SPOT-VGT data (1982–2004). Remote Sensing of Environment, 101, 52−62. Edwards, M., & Richardson, A. J. (2004). Impact of climate change on marine pelagic phenology and trophic mismatch. Nature, 430, 881−884. ESRI (Environmental Systems Research Institute) (2006). ArcMap 9.2. Redlands, California: ESRI. Fagre, D. B., Peterson, D. L., & Hessel, A. E. (2003). Taking the pulse of mountains: Ecosystems response to climate variability. Climate Change, 59, 263−282. Food and Agriculture Organization (FAO) of the United Nations, Economic and Social Department, 2000. Friedl, M. A., & Brodley, C. E. (1997). Decision tree classification of land cover from remotely sensed data. Remote Sensing of Environment, 61, 399−409. Friedl, M. A., McIver, D. K., Hodges, J. C. F., Zhang, X. Y., Muchoney, D., Strahler, A. H., et al. (2002). Global land cover mapping from MODIS: Algorithms and early results. Remote Sensing of Environment, 83, 287−302. Goetz, S. J., Bunn, A. G., Fiske, G., & Houghton, R. A. (2005). Satellite observed photosynthetic trends across boreal North America associated with climate and fire disturbance. Proceedings of the National Academy of Science, 102(38), 13521−13525. Goetz, S. J., & Prince, S. D. (1996). Remote sensing of net primary production in boreal forest stands. Agriculture and Forest Meteorology, 78(3), 149−179. Gu, L., Hanson, P. J., Mac Post, W., Kaiser, D. P., Yang, B., Nemani, R., et al. (2008). The 2007 eastern US spring freeze: Increased cold damage in a warming world. BioSciences, 58, 253−262. Gu, L., Meyers, T., Pallardy, S. G., Hanson, P. J., Yang, B., Heuer, M., et al. (2007). Influences of biomass heat and biochemical energy storages on the land surface fluxes and radiative temperature. Journal of Geophysical Research, 112, d02107. doi:10.1029/2006JD007425. Hardisky, M. A., Klemas, V., & Smart, R. M. (1983). The influence of soil, salinity, growth form, and leaf moisture on the spectral radiance of Spartina alterniflora canopies. Photogrammertric Engineering and Remote Sensing, 49, 77−83. Hopkins, A. D. (1918). The bioclimatic law. Monthly Weather Review(Suppl. 9). Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83, 195−213. Huete, A., Didan, K., Shimabukuro, Y. E., Ratana, P., Saleska, S. R., Hutyra, L. R., et al. (2006). Amazon rainforests green-up with sunlight in dry season. Geophysical Research Letters, 33. doi:10.1029/2005GL025583. Hunt, E. R., & Rock, B. N. (1989). Detection of changes in leaf water content using near- and middle-infrared reflectances. Remote Sensing of Environment, 30, 43−54. Inouye, D. W. (2008). Effects of climate change on phenology, frost damage, and floral abundance of mountain wildflowers. Ecology, 89, 353−362.

Inouye, D. W., & Wielgolaski, F. E. (2003). Phenology of high-altitude climates. In M. D. Schwartz (Ed.), Phenology: An integrative environmental science (pp. 195−214). Netherlands: Kluwer Academic publishers. Ives, J. D., Messerli, B., & Spiess, E. (1997). Mountains of the world — A global perspective. In B. Messerli, & J. D. Ives (Eds.), Mountains of the world: A global priority (pp. 313−335). Oxford: Parthenon Publication Group. Jiang, Z., Huete, A., Didan, K., & Miura, T. (2008). Development of a two-band enhanced vegetation index without a blue band. Remote Sensing of Environment, 112, 3833−3845. Körner, C. (2005). The green cover of mountains in a changing environment. In U. M. Huber (Ed.), Global change in mountain regions (pp. 367−375). Netherlands: Springer Press. Lieth, H. (1974). Phenology and seasonality modeling. New York: Springer-Verlag. Lucht, W. (1998). Expected retrieval accuracies of bi-directional reflectance and albedo from EOS-MODIS and MISR angular sampling. Journal of Geophysical Research, 103, 8763−8778. Lucht, W., Hyman, A. H., Strahler, A. H., Barnsley, M. J., Hobson, P., & Muller, J. -P. (2000). A comparison of satellite-derived spectral albedos to ground-based broadband albedo measurements modelled to satellite spatial scale for a semi-desert landscape. Remote Sensing of Environment, 74, 85−98. Lucht, W., & Lewis, P. (2000). Theoretical noise sensitivity of BRDF and albedo retrieval from the EOS-MODIS and MISR sensors with respect to angular sampling. International Journal of Remote Sensing, 21, 81−98. Lucht, W., & Roujean, J. -L. (2000). Considerations in the parametric modeling of BRDF and Albedo from multiangular satellite sensor observations. Remote Sensing Review, 18, 343−379. Lucht, W., Schaaf, C. B., & Strahler, A. H. (2000). An algorithm for the retrieval of albedo from space using semiempirical BRDF models. IEEE Transactions in Geosciences, Remote Sensing, 38, 977−998. Malingreau, J. P. (1986). Global vegetation dynamics: Satellite observations over Asia. International Journal of Remote Sensing, 7, 1121−1146. Mark, A. F. (1970). Floral initiation and development in New Zealand alpine plants. New Zealand Journal of Botany, 8, 67−75. Meehl, G. A., Karl, T., Easterling, D. R., Changnon, S., Pielke, R., Changnon, D., et al. (2000). An introduction to trends in extreme weather and climate events: Observations, socioeconomic impacts, terrestrial ecological impacts, and model predictions. Bulletin of the American Meteorological Society, 81, 413−416. Menzel, A., & Fabian, P. (1999). Growing season extended in Europe. Nature, 397, 659. Morisette, J. T., Richardson, A. D., Knapp, A. K., Fisher, J. I., Graham, E., Abatzoglou, J., et al. (2009). Unlocking the rhythm of the seasons in the face of global change: Challenges and opportunities for phenological research in the 21st century. Frontiers in Ecology and the Environment, 7, 253−260. Moulin, S., Kergoat, L., & Viovy, N. (1997). Global-scale assessment of vegetation phenology using NOAA/AVHRR satellite measurements. Journal of Climate, 10, 1154−1170. Muchoney, D., Borak, J., Chi, H., Friedl, M. A., Hodges, J., Morrow, N., et al. (2000). Application of the MODIS global supervised classification model to vegetation and land cover mapping of Central America. International Journal of Remote Sensing, 21, 1115−1138. Myneni, R. B., Keeling, C. D., Tucker, C. J., Asrar, G., & Nemani, R. R. (1997). Increased plant growth in the northern high latitudes from 1981 to 1991. Nature, 386, 698−702. National Climate Archives, New Brunswick, Canada. 2007. http://www.climate. weatheroffice.gc.ca/climateData/dailydata_e.html?StationID=10223 Pachauri, R. K., & Reisinger, A. (Eds.). (2007). Climate change 2007 synthesis report (pp. 104). Geneva, Switzerland: IPCC. Parmesan, C., & Yohe, G. (2003). A globally coherent fingerprint of climate change across natural systems. Nature, 421, 37−42. Peñuelas, J., & Filella, I. (2001). Responses to a warming world. Science, 294, 793−795. Post, E. (2003). Timing of reproduction in large mammals. In M. D. Schwartz (Ed.), Phenology: An integrative environmental science (pp. 195−214). Netherlands: Kluwer Academic Publishers. Price, L. W. (1981). Mountains and man: A study of processes and environment. Berkeley, London: University of California Press. Price, M. F. (1990). Mountains — globally important ecosystems. Unasylva, 49. PRISM Climate Group. 2010. Oregon State University, http://www.prismclimate.org Ratcliffe, M. J., & Turkington, R. (1989). Comparative phenology of some alpine vascular plant species on Lakeview Mountain, southern British Columbia. Canadian Field Naturalist, 103, 348−352. Reed, B. C., Brown, J. F., VanderZee, D., Loveland, T. R., Merchant, J. W., & Ohlen, D. O. (1994). Measuring phenological variability from satellite imagery. Journal of Vegetation Science, 5, 703−714. Rich, P. M., Dubayah, R., Hetrick, W. A., & Saving, S. C. (1994). Using Viewshed models to calculate intercepted solar radiation: Applications in ecology. American Society for Photogrammetry and Remote Sensing Technical Papers (pp. 524−529). Roy, D. P., Ju, J., Kline, K., Scaramuzza, P. L., Kovalskyy, V., Hansen, M., et al. (2009). Webenabled Landsat Data (WELD) preliminary results: Landsat ETM + composited mosaics of the conterminous United States. Remote Sensing of Environment, 114, 35−49. Schwartz, M. D., & Reiter, B. E. (2000). Changes in North American spring. International Journal of Climatology, 20, 929−932. Schwartz, M. D., Ahas, R., & Anto, A. (2006). Onset of spring starting earlier across the Northern Hemisphere. Global Change Biology, 12, 343−351. Sphen, E. M., & Korner, C. (2002). A global assessment of mountain biodiversity and its function. In C. Korner, & E. M. Sphen (Eds.), Mountain biodiversity: A global assessment (pp. 393−400). Oxford: Parthenon Publication Group.

A.D. Hudson Dunn, K.M. de Beurs / Remote Sensing of Environment 115 (2011) 1220–1233 Tommervik, H. (2004). Vegetation changes in the Nordic mountain birch forest: The influence of grazing and climate change. Arctic Antarctic and Alpine Research, 36, 323−332. Tucker, C. J. (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment, 8, 127−150. Tucker, C. J., Slayback, D. A., Pinzon, J. E., Los, S. O., Myneni, R. B., & Taylor, M. G. (2001). Higher northern latitude normalized difference vegetation index and growing season trends from 1982 to 1999. International Journal of Biometeorology, 45, 184−190. Tucker, C. J., Townshend, J. R. G., & Goff, T. E. (1985). African land-cover classification using satellite data. Science, 227(4685), 369−375. United States Department of Agriculture, National Resources Conservation Service (http://www.wcc.nrcs.usda.gov/snotel/) Wagner, J., & Reicheggar, B. (1997). Phenology and seed development of the alpine sedges Carex curvula and Carex firma in response to contrasting topoclimates. Arctic and Alpine Research, 29, 291−299. Walker, M. D., Ingersoll, R. D., & Webber, P. J. (1995). Effects of inter-annual climate variation on phenology and growth of two alpine forbes. Ecology, 76, 1067−1083. Walsh, S. J. (2004). An assessment of snow avalanche paths and forest dynamics using IKONOS satellite data. Geocarto International, 19, 85−93. Weiss, D. J., & Walsh, S. J. (2009). Remote sensing of mountain environments. Geography Compass, 3(1), 1−21. White, M. A., de Beurs, K. M., Kamel, D., Inouye, D. W., Richardson, A. D., Jensen, O. P., et al. (2009). Intercomparison, interpretation, and assessment of spring phenology

1233

in North America estimated from remote sensing for 1982 to 2006. Global Change Biology, 15, 2335−2359. White, M. A., Hauffman, F., Hargrove, W. W., & Nemani, R. R. (2005). A global framework for monitoring phenological responses to climate change. Geophysical Research Letters, 32(L0), 4705. doi:10.1029/2004GL021961. White, M. A., Thornton, P. E., & Running, S. W. (1997). A continental phenology model for monitoring vegetation responses to interannual climatic variability. Global Biogeochemical Cycles, 11, 217−234. Xiao, X., Boles, S., Frolking, S., Sala, W., Moore, B., III, & Li, C. (2002). Observation of flooding and rice transplanting of paddy rice fields at the site to landscape scales in China using VEGETATION sensor data. International Journal of Remote Sensing, 23 (15), 3009−3022. Xiao, X., Boles, S., Liu, J., Zhuang, D., & Liu, M. (2002). Characterization of forest types in Northeastern China, using multi-temporal SPOT-4 VEGETATION sensor data. Remote Sensing of Environment, 82, 335−348. Zhang, X. Y., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C., Gao, F., et al. (2003). Monitoring vegetation phenology using MODIS. Remote Sensing of Environment, 84, 471−475. Zhang, X. Y., Friedl, M. A., Schaaf, C. B., Strahler, A. H., & Schneider, A. (2004). The footprint of urban climates on vegetation phenology. Geophysical Research Letters, 31(L1), 2209. doi:10.1029/2004GL020137.