Spatial patterns of land surface phenology relative to monthly climate variations: US Great Plains

Abstract

We extracted and mapped six land surface phenological metrics including: (1) the peak normalized difference vegetation index (NDVI), (2) peak date, (3) start of season (SOS), (4) end of season (EOS), (5) length of growing season (LOS), and (6) cumulative NDVI from 2000 to 2009 using Moderate-Resolution Imaging Spectroradiometer (MODIS) images covering the United States (US) Great Plains. Their patterns relative to monthly precipitation, maximum temperature, minimum temperature, and dew points were analyzed using multiple linear regression, stepwise selection, and geographically weighted regression (GWR) analysis. Both peak NDVI and cumulative NDVI had similar spatial patterns. Their values decreased along an east to west gradient. Peak date and SOS also showed compatible patterns. The southeastern Great Plains had the earliest SOS, peak date, and the longest LOS, given its warmer temperatures and greater precipitation. Dew points in March and October as well as the maximum temperature in April highly influenced the SOS, while dew point in August was found more influential for EOS and LOS. Precipitation in March and September also affected the total cumulative NDVI. The GWR models performed better than the OLS because the GWR utilized the spatial relationships between the different variables resulting from local level processes. The regression models predicted peak NDVI and cumulative NDVI better than the other phenological indices.

Keywords:

• land surface phenology
• the US Great Plains
• climate variations
• MODIS
• geographically weighted regression

1. Introduction

Plant phenology can be defined as the study of the timing of plant’s recurring cycles as influenced by the environment and weather conditions (Cleland et al. 2007). Global climate change could significantly alter plant phenology by influencing the timing of growth development and driving annual carbon uptake in terrestrial ecosystems (Partanen, Koski, and Hänninen 1998; Cleland et al. 2007; Sanz-Perez, Castro-Diez, and Valladares 2009). Hansen et al. (2006) found that global temperature increased 0.7°C from the late nineteenth century to 2000, and continued warming in the twenty-first century remains persistent with an increase of 0.2°C in the past decade. In this context, phenology shift has been identified at multiple (species, ecosystem and global) scales (Cleland et al. 2007). For example, Kunkel et al. (2004) found an increase of approximately 2 weeks in frost-free season length in the twentieth century, which conforms to the findings of Schwartz, Ahas, and Aasa (2006) that the first leaf dates and spring last frost dates were 1.2 and 1.5 days earlier per decade, respectively, in the temperate Northern Hemisphere from 1955 to 2002. Miller, Mitchell, and Lopez (2005) reported a similar finding relative to the start of the growing season in the US Corn Belt. Consequently, the acceleration of spring has resulted in a longer growing season (Robeson 2004). In Central Europe, spring time frost events advanced faster than phenological phases from 1951 to 1997 (Scheifinger et al. 2003). These phenoclimatic measures have important indications for agricultural practices and parameterizing global climate models (Cleland et al. 2007). Phenology also influences seasonal fluxes of water, energy, and carbon between the land surface and atmosphere (Ganguly et al. 2010). It is important to characterize the global phenology because of strong interactions between climate and phenological processes in almost all ecosystems exhibiting seasonal weather variability (Parmesan and Yohe 2003; Hudson-Dunn and de Beurs 2011).

Different from field species-level phenological observations, remote sensing provides a unique way to monitor land surface phenology (LSP) over large scales, varying land covers and at regular intervals. LSP studies the spatiotemporal patterns of the vegetated land surface as observed by synoptic sensors at spatial resolutions and extents relevant to meteorological processes (de Beurs and Henebry 2005). Information related to LSP is widely used as a diagnostic measure of ecosystem response to global change (Ganguly et al. 2010). In the past three decades, various remote-sensing methods using daily vegetation indices (VI) collected by the Advanced Very High Resolution Radiometer (AVHRR), the Moderate Resolution Spectroradiometers (MODIS), or SPOT-VEGETATION sensors have been developed for monitoring seasonal vegetation changes (Justice  et al. 1985; White, Thornton, and Running 1997; Xiao et al. 2002; Kang et al. 2003; Zhang et al. 2003; Delbart et al. 2006).

To gain better understanding of the effects of climate change on plant phenology, the relationships between normalized difference vegetation index (NDVI) and climatic variables have been widely investigated. The monthly timing of NDVI extremes was found closely associated with seasonal patterns of maximum and minimum temperature and rainfall (Potter and Brooks 1998). Although temperature was found to be the primary climatic factor associated with greening in the northern high latitudes, precipitation proved a strong correlate of vegetation growth in fragmented or arid regions (Nicholson, Davenport, and Malo 1990; Nicholson and Farrar 1994; Eklundh 1998; Du Plessis 1999; Kawabata, Ichii, and Yamaguchi 2001; Wang, Rich, and Price 2003; Xiao and Moody 2005). In most temperate regions, vegetation dynamics are affected both by precipitation and temperature (Schultz and Halpert 1993). In most of these previous studies, multiple linear regression models related to NDVI and climatic factors were constructed. These models are based on the global technique of ordinary least squared (OLS) regression analysis with the assumption of spatial stationarity. However, the use of this global technique over relatively large areas is constrained by local spatial variability of the observed relationships.

To take these local constraints into consideration, spatial statistical approaches, such as geographically weighted regression (GWR), become more appropriate (Maselli 2002). Foody (2003) revealed the potential use of GWR to the remote-sensing community, showing that GWR may present a more effective way of analyzing spatially varying relationships between the NDVI and rainfall over North Africa and the Middle East. Using OLS and GWR, Foody and Dash (2010) explored the relationship between the Envisat Medium-Resolution Imaging Spectrometer (MERIS) Terrestrial Chlorophyll Index (MTCI) and C3 grassland composition (C₃%) of the Great Plains. Their results indicated that a relationship between MTCI data and grassland composition may be formulated in the State of South Dakota with R² = 0.6, but the transferability of the relationship to other states was limited by its significant nonstationarity in a large extent. Roy and Yuan (2007) examined the relationship between summer NDVI and summer precipitation and temperature in Minnesota from 1990 to 1997 using GWR. They found overall relationships between NDVI and the two climatic factors were not globally stationary, but the spatial patterns of the correlation were consistent with different types of land uses and prevailing vegetation across the state. They also indicated that the GWR model performance (average R² = 0.67) would be more significant if environmental factors such as soil moisture could be added in future studies. Finally, a GWR analysis was used successfully to document localized vegetation recovery in the Black Hills of South Dakota following a major forest fire that occurred in 2000 (Mitchell and Yuan 2010).

Drawing upon these previous studies, this research aims to identify the significant climate factors for phenological change and their local variations in the Great Plains using OLS and GWR models. The specific objectives are twofold: first, six land surface phenological metrics – the peak NDVI, peak date, start of season (SOS), end of season (EOS), length of growing season (LOS), and cumulative NDVI – were extracted from the 2000–2009 NDVI time series. Their spatial and temporal patterns were examined thereafter. Second, the relationships between these phenological metrics and monthly climate variables (precipitation, maximum temperature, minimum temperature, and dew points) were analyzed using the OLS and GWR regression models. The results help provide a better understanding of ecosystem dynamics in this region within the context of global climate change.

2. Study area

The US Great Plains lie west of the Mississippi River and east of the Rocky Mountains in the United States and possess a mid-latitude steppe climate characterized by a water balance deficit with precipitation and soil moisture supplying between 50% to less than 100% of potential evapotranspiration. The climate also has large seasonal temperature ranges due to its low specific heat associated with a continental location (Oliver and Hidore 2002). Covering 10 states – Montana, Wyoming, Colorado, New Mexico, North Dakota, South Dakota, Nebraska, Kansas, Oklahoma, and Texas – the core area of Great Plains extends about 800 km (500 miles) east to west and 3200 km (2000 miles) north to south (Thornbury 1965; Christopherson 1997).

As the primary grassland of North America, the Great Plains contain approximately 162 million ha (625,000 sq. miles) of prairies, much of which has been altered by Euro–American agriculture. The main prairie bodies within the Great Plains are: (1) the tall grass prairie extending from Canada and Minnesota south to Texas, (2) the mixed grass prairie from Canada and eastern North Dakota south to Texas, and (3) the short grass plains reaching from Montana to the Llano Estacado of west Texas and eastern New Mexico. Although parts of the short grass prairie have been modified by irrigation agriculture (as in the vicinity of Lubbock, Texas) and elsewhere the dryland farming of wheat is pervasive; a majority of short grass prairie presently remains rangeland grazed by cattle, a situation mimicking the native bison herds. The mixed grass prairie includes a combination of crop and range land with more limited stands of native vegetation as exhibited at Badlands National Park, South Dakota. Since 1830s, the majority of tall grass prairie underlain by fertile and highly productive mollisol soils has been converted to cropland. The declines in area of tall grass prairie (82–99%) exceed those reported for any other major ecosystems in North America (Sampson and Knopf 2012).

Aside from the variations of floristic regions in these prairies, a noticeable latitudinal variation in plant functional types (PFTs) exists in the Great Plains. Cold-season C3 and warm-season C4, which differ from their photosynthesis pathways, are the two most widely distributed PFTs present in the Great Plains. Relying on growing-season temperature, grasslands in the southern Great Plains are dominated by C4 species, while C3 species become increasingly dominant northwards (Tieszen et al. 1997; Foody and Dash 2010). Wang et al. (2013) investigated the relative abundance of C3 and C4 PFTs throughout the Great Plains. C3 and C4 photosynthesis are relevant to global change studies since they respond very differently to changes in climatic conditions. C4 favors higher growing-season temperature and lower atmospheric CO₂, while the reverse is true for the C3. With its great latitudinal variation and the different grassland compositions, the Great Plains provides an ideal case study site for analyzing local phenological patterns and their dependences to climatic variations.

3. Methodology

3.1. Data sets

The core study area is covered in four MODIS tiles (Figure 1).  For each tile, a total of 46 MOD09A1 surface reflectance composites (8-day, 500-m) were downloaded from 2000 to 2009 from the Land Processes Distributed Active Archive Center (http://lpdaac.usgs.gov). Pouliot et al. (2011) found composite periods from 7 to 11 days produced the best results for some LSP metrics detection. The NDVI time series each year were then extracted from these image products using the band 1 (red) and band 2 (near-infrared) of the MOD09A1 products. All pixels with NDVI < 0.1 were treated as nonvegetation covers because they were assumed photosynthetically inactive in land surface phenology. Because the MOD09A1 products have already been georectified and the images acquired in different dates matched well, additional geometric rectification of MODIS time series was not performed.

Figure 1. The US Great Plains in the four tiles of July 2007 MODIS images (bands 2, 4, and 3 false color composite).

Monthly climate data (precipitation, max. temperature, min. temperature, and dew point) for the 10-year duration were extracted from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) data. The PRISM is a knowledge-based system that uses point meteorological measures at ground stations to produce continuous, digital grid estimates of climatic parameters (http://prism.oregonstate.edu/). A 2007–2008 Cropland Data Layer (CDL) land cover classification map acquired from the US Department of Agriculture (USDA) was modified by combining some of the crop types and used as a visual reference map for our results (Figure 2).

3.2. Phenological metrics extraction

The procedure of phenological metrics extraction includes three general steps: data filtering, temporal smoothing, and derivation of metrics. In particular, the extracted eight-day NDVI time series has been processed with the maximum-value-composite (MVC) technique to remove clouds, haze, and snow pixels. The MVC method assumes the presence of clouds, haze, and snow in a pixel will reduce the NDVI value. Hence, each NDVI value is examined on a pixel-by-pixel basis. The highest value of the eight days is retained for each pixel. However, the 46-scene NDVI time series from each year are still jagged, because of cloud residuals and varying aerosol conditions in different dates. Hence, a five-point median filter was applied to remove spikes and troughs in the trajectories (Wang et al. 2011). Next, a Savitzky–Golay filtering method (Savitzky and Golay 1964) was used to smooth out minor atmospheric effects in the time series curves (Wang et al. 2011). The Savitzky–Golay smoothing filter performs a local polynomial regression on a series of values to determine the smoothed value for each point. It preserves features of the distribution such as relative maxima, minima, and width, which are usually ‘flattened’ by other adjacent averaging techniques.

With the 46-scene NDVI time series, a thresholding approach was applied to extract the six phenology metrics at each pixel, which included (1) start of season (SOS): the date (day-of-year, or DOY) when NDVI increases to 20% of the amplitude; (2) end of season (EOS): the date when NDVI decreases to 20% of the amplitude; (3) LOS: the number of dates between SOS and EOS; (4) peak NDVI: the maximum NDVI value in a growing season; (5) peak date: the date that peak NDVI falls on; (6) cumulative NDVI: the integral of NDVI from SOS to EOS, which approximates biomass accumulation of each pixel. For SOS and EOS, a threshold of 20% was selected because it was commonly applied in studies of grassland phenology (Zhang et al. 2003; Wang et al. 2011).

3.3. Multivariate regressions and GWR analysis

With a reference map (Figure 2), a total of 2000 vegetation samples were randomly selected across the study area. From 2000 to 2009, phenological metrics and climate data were extracted. Their relationships in each year were first examined using the global multiple linear regression OLS method. Then, the GWR analysis was performed to explore the local variations of spatial patterns of the relationships. Equation (1) describes the OLS model:

(1)

where the dependent variable was one of the six phenological metrics (PM); the 50 independent variables were monthly maximum temperature (T_Max), minimum temperature (T_Min), precipitation (P), and dew point temperature (D) from January to December, and the mean maximum (T_{Max_mean}) and minimum temperatures (T_{Min_mean}) of the modeled year. This model assumed the spatial stationarity in the relationship among the 50 variables under study, where β₀ is the intercept and β₁ to β₅₀ are the slopes; and ε is the residual.

A stepwise regression with a bidirectional elimination approach was applied to determine the final predictive variables from the 50-candidate independent variables in Equation (1). The identified independent variables were assumed to have statistically significant correlation with the dependent variable and were then used in the multivariate GWR analysis. In comparison to the OLS multiple regression, the base equation used in the multivariate GWR was described in Equation (2):

(2)

where (u,v) are the geographic coordinates of a particular point in the study area. In GWR, the model parameters (β) are assumed to vary locally, which are calculated as a weighted contribution in accordance to its spatial proximity to the specific location being considered. Therefore, the derived parameters can be mapped to identify the nature of their variation in space to bring out the spatially varying relationships among the different variables. In this sense, GWR is an effective method to examine the presence of spatial nonstationarity in relationships (Fotheringham, Charlton, and Brundson 1998). Because the 2000 sample points were randomly and densely distributed across the study area, a fixed kernel method with Gaussian function was conducted. The weight of the jth data point at the ith regression point is calculated by following equation:

(3)

where d_ij is the Euclidean distance between the regression point i and the data point j, and b is the bandwidth. The weights of data points decrease as the distance from the regression point increases. The bandwidth (b) was selected according to the Akaike Information Criteria (AIC) minimization as defined by Equation (4) (Fotheringham, Brunsdon, and Charlton 2002):

(4)

where n is the sample size (n = 2000 in this study), is the estimated standard deviation of the error term, and tr(S) is the trace of the hat matrix S, which projects (the modeled phenological index) onto (the observed phenological index) by .

Figure 2. Land covers in the US Great Plains (modified from the 2007/08 CDL map).

Finally, an ANOVA tested the null hypothesis that no improvement exists between the GWR and the OLS model. Furthermore, to illustrate the overall spatial patterns across the study area during the 10-year span, the local coefficients of determination (R²) of the GWR analysis were interpolated into continuous surfaces using an Ordinary Kriging interpolation method, which assumes an unknown constant trend.

4. Results and discussion

4.1. Climate patterns

Monthly climate demonstrated dramatic temperature changes from winter to summer (Figure 3). The mean minimum temperature values of the study area ranged from –10°C in December–January to 17.5°C in July, with standard deviations varying from 3.1°C to 7.2°C during the 10-year span (Figure 3(a)). For the monthly maximum temperature, the ranges were 3°C to 33°C for the mean values and 1.8°C to 8.6°C for the standard deviations (Figure 3(b)). The standard deviation of the maximum temperature showed apparent seasonality in the 10-year period and revealed much higher variation in the winters and more spatially stable (lower variation) in growing seasons. These findings indicate greater air mass variability in the winter as the region or portions thereof can either be impacted by: (1) highly modified maritime polar air masses associated with a zonal jet (west to east), or (2) blasts of cold, dense, and dry continental polar air in the context of a more meridional jet stream (west to east flow with strong north–south sinuosity).

Similar to that of the minimum temperature, the dew point standard deviation did not reveal clear seasonality in the 2000–2009 timespan. On the other hand, erratic monthly and yearly precipitation changes with large variance are displayed in Figure 3(d), which indicate strong variations in synoptic conditions associated with enhanced precipitation or to the contrary drought across the study area. For example, much of the Great Plains contains a bimodal set of precipitation controls. Wintertime totals are lower and associated with frontal activity (cyclogenesis) that can either be enhanced or absent depending on jet stream configurations. Similarly, a late spring or summer time regime possesses variability based on the extent and intensity of maritime tropical air penetration emanating from the Gulf of Mexico, which accentuates frontal activity or results in convectional air mass thunderstorms. The Dust Bowl years of the 1930s were for instance noted for failures of both seasonal precipitation regimes on a repetitive basis and resulted in soil moisture deficits of 34 inches (0.8 m) in western Kansas (Worster 2004).

Figure 3. Average and standard deviations of monthly weather conditions of the Great Plains from 2000 to 2009.

From the spatially averaged temperature imagery, monthly maximum temperature decreases from south to north, as expected, following trends in net absorbed solar radiation. Monthly minimum temperature and dew point demonstrated a southeast to northwest gradient (Figures 4–6), which again follows a trajectory of maritime tropical penetration from the Gulf of Mexico, the primary source of moisture for much of the US east of the Rockies. Indicative of low specific heat and continentality, large temperature ranges (>30°C) were shown in almost all the monthly maximum and minimum temperatures. Spatial patterns of significant variability in precipitation were revealed in Figure 7. As noted, much of the Great Plains lies in the semi-arid mid-latitude steppe climate, where annual precipitation is greater than half but less than 100% of the annual potential evapotranspiration. Interior continental and leeward orographic position relative to the Sierra-Cascade and Rocky Mountains in part explain this precipitation deficiency. The Rockies also influence jet stream dynamics, which as noted above can either accentuate or greatly minimize precipitation events. Location relative to the extent, intensity and duration of maritime tropical air from the Gulf of Mexico stands as yet another factor that can either aridify or enhance water availability (Oliver and Hidore 2002). Indeed, the 10-year mean monthly precipitation ranged from 0 mm to 217 mm with a positive and sharp precipitation gradient decreasing from the east toward the west in the growing season months (Figure 7). For example, the highest totals of annual precipitation (more than 1300 mm) were located along or close to 97°W in Texas, Oklahoma, and southern Kansas, in areas where maritime tropical air masses have a higher probability of occurrence.

Figure 4. Ten-year mean monthly maximum temperature (°C) in the US Great Plains.

Figure 5. Ten-year mean monthly minimum temperature (°C) in the US Great Plains.

Figure 6. Ten-year mean monthly dew point (°C) in the US Great Plains.

Figure 7. Ten-year mean monthly precipitation (mm) in the US Great Plains.

4.2. Phenological metrics patterns

The images of six phenological indices for the Great Plains were derived based on the smoothed NDVI time series. The average spatial patterns of the 10 years were then calculated and displayed in Figure 8. We can see the mean peak NDVI (Figure 8(a)) and the mean cumulative NDVI (Figure 8(f)) displayed similar spatial patterns. Their values decreased in an east–west direction. Eastern tall grass prairie had the highest peak and cumulative NDVIs associated with greater and more reliable precipitation. Peak date and the SOS (Figure 8(b) and 8(c)) exhibited a comparable pattern. The ending point of growing season (EOS) proved more temperature-driven across the whole of the study area. Therefore, a more homogeneous spatial pattern of EOS was produced (Figure 8(d)).

Figure 8. Ten-year mean phenological metrics extracted from multitemporal MODIS data. (a) Peak NDVI; (b) Peak date (DOY); (c) SOS (DOY); (d) EOS (DOY); (e) LOS (Days); and (f) Cumulative NDVI.

The southeastern Great Plains possessed the earliest vegetation peak date and the SOS, as well as the longest LOS, as determined by their warmer temperature and relatively wetter conditions. Conversely, the southwestern portion of the Great Plains (eastern New Mexico and northwestern Texas) has the latest SOS. For example, the southern plains of Texas (Panhandle region) peak in precipitation during June or late spring with a noticeable reduction in July and August. The areas near Amarillo and Lubbock, Texas, fit this mode of precipitation. As such, the warm season short grasses will green-up during the late spring and early summer.

Meanwhile, in eastern New Mexico, precipitation peaks during August and September due to the Arizona Monsoon, which results from a thermal low developing over southern Arizona and New Mexico, and draws moisture inland from the Gulf of California and Gulf of Mexico, respectively (Rohi and Vega 2008). Both Tucumcari and Hobbs, New Mexico fit this pattern and are located on the most downwind side of the Arizona Monsoon. Consequently, warm season short grasses green-up in the late summer (late August through September). As Figure 9 demonstrates, large SOS interannual variations occurred in the southwest corner of the Plains from 2000 to 2009. This SOS situation again points to the extreme variability of the region’s precipitation regime, a factor that is supplemented by (1) the reliability of long hot summers, and (2) a change in climate towards a shorter nongrowing season, which enhances early spring evapotranspiration and exasperates water balance deficits in an already marginal area bordering true desert climates to the south.

Figure 9. Normalized difference map of SOS from 2000 to 2009.

Figure 10 highlights the Cumulative NDVI relative to the 10-year mean also indicated more interannual variability occurred in the tall grass prairie section of the Great Plains. Although this finding was not as strong as the aforementioned SOS interannual variability, it follows the general trend associated with a warmer atmosphere, namely, a greater propensity for heavier rainfall events separated by longer dry spells.

Figure 10. Normalized difference map of cumulative NDVI from 2000 to 2009.

4.3. Effects of climate factors on phenological metrics

Spatial distributions of these phenological indices in the Great Plains revealed climate dependences of vegetation growth in this large region. In order to determine the climate variables that most significantly impact vegetation phenology, the stepwise regression results for all the OLS models are summarized in Table 1. These variables were selected based on a threshold of 0.05 significance level for at least 5 of the 10 years. Dew point played the most significant role among all the climate variables except for the EOS. For the Peak NDVI, dew point was significant in almost all months across a year (March, June, December, February, November, August, and September; sorted by its significance in descending order), while precipitation in May and February were significant. Dew point temperatures in February, December, June, and November and precipitation in February and May significantly determined the peak NDVI values. Cumulative NDVI was significantly influenced by precipitation in early spring (March) and fall (September and October), which effectively served to load the soil with adequate moisture. For Peak date, spring (March and May) and October dew points were the most significant factors. For SOS, January and October’s dew points were the most important variables. The dew points suggest more water vapor present, itself a highly potent greenhouse gas that would limit nocturnal cooling. From a practical standpoint, it is difficult to reason why January would test significant given its standing as being in the midst of non-growing season where soils are mostly frozen especially in the northern plains. Perhaps this situation is a statistical quirk implying the need to use the GWR. Meanwhile, the EOS proved highly affected by November’s precipitation. We postulate that cloud cover associated with enhanced precipitation and frontal activity limited nocturnal temperatures. Finally, the LOS was greatly determined by the dew point in March. Again, the presence of a greenhouse gas, in this case water vapor, would limit nocturnal long wave radiation losses and thereby favor warmer nighttime temperatures.

Table 1. Most significant predictive variables identified in the OLS models during the 10-year timespan. The variable is significant at the 0.05 significance level for at least 5 years out of the 10 year-span.

Dependent variable	Most significant independent variables identified in the OLS models from 2000 to 2009
Peak NDVI	dew3, dew6, dew12, pre5, dew2, dew11, dew8, dew9, pre2
Peak date	dew3, dew10, dew5
SOS	dew1, dew10
EOS	pre11
LOS	dew3
Acc. NDVI	dew2, dew12, pre3, pre10, dew6, dew11, pre9

The results of the OLS and GWR analyses are summarized in Tables 2 and 3. Table 2 lists the coefficient of determination (R²) values of the OLS and their plain-level averages in the GWR analyses. The highest correlations with R² greater than 0.5 for all the years both for OLS and GWR analyses were found for Peak NDVI and Cumulative NDVI when used as the dependent variables. Comparatively, both models explained the lowest amount of variance for EOS except in 2000, suggesting a more stable EOS. This finding conforms to the Miller, Mitchell, and Lopez (2005) study of growing season length in the US Corn Belt that was based on recorded temperatures rather than satellite data. Our OLS models from 2000 to 2009, by using monthly precipitation, temperatures, as well as dew point values to model six phenological indices, performed better than the seasonal based OLS model of Roy and Yuan (2007). Furthermore, the GWR model always reached higher R² than the OLS for all phenological indices, indicating that the GWR better explained climate dependencies of vegetation phenology by considering both global trends and local variations.

Table 2. Comparison of results (R²) from OLS and GWR models.

Dependent variable	2000		2001		2002		2003		2004		2005		2006		2007		2008		2009
	OLS	GWR	OLS	GWR	OLS	GWR	OLS	GWR	OLS	GWR	OLS	GWR	OLS	GWR	OLS	GWR	OLS	GWR	OLS	GWR
Peak NDVI	0.63	0.73	0.65	0.75	0.65	0.73	0.59	0.71	0.63	0.74	0.62	0.71	0.57	0.68	0.59	0.70	0.62	0.72	0.64	0.72
Peak date	0.28	0.36	0.31	0.36	0.21	0.30	0.28	0.38	0.20	0.32	0.20	0.28	0.31	0.39	0.29	0.41	0.32	0.44	0.27	0.32
SOS	0.22	0.30	0.27	0.38	0.33	0.37	0.22	0.29	0.23	0.28	0.27	0.35	0.28	0.39	0.34	0.38	0.38	0.51	0.25	0.33
EOS	0.15	0.32	0.16	0.40	0.08	0.17	0.14	0.28	0.13	0.26	0.11	0.21	0.15	0.26	0.27	0.43	0.13	0.30	0.09	0.25
LOS	0.09	0.28	0.20	0.42	0.21	0.29	0.18	0.31	0.21	0.35	0.15	0.23	0.17	0.30	0.30	0.44	0.21	0.38	0.11	0.26
Acc. NDVI	0.51	0.67	0.48	0.67	0.53	0.64	0.51	0.66	0.50	0.66	0.58	0.69	0.54	0.66	0.58	0.70	0.58	0.69	0.58	0.68

The ANOVA analysis results and the F-values shown in Table 3 confirm that the GWR models are a significant improvement on the global OLS models with all dependent variables, indicating apparent spatial nonstationarity exists in the relationships between the phenological metrics and the monthly weather conditions. This result agrees with the earlier study by Foody and Dash (2010) in the Great Plains.

Table 3. F-values of the ANOVA analysis.

Dependent variable	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009
Peak NDVI	4.44	3.87	3.94	3.83	3.61	3.38	2.69	3.73	3.78	3.74
Peak date	2.75	3.98	3.17	3.52	2.88	3.08	3.90	3.26	4.26	3.46
SOS	3.52	3.89	2.79	2.88	3.11	3.25	3.62	3.03	2.93	3.22
EOS	3.18	4.18	2.36	3.42	3.18	3.20	4.04	3.20	3.99	3.17
LOS	3.44	4.31	2.47	3.17	3.59	3.45	3.70	3.11	3.47	3.05
Acc. NDVI	4.11	4.74	4.01	4.09	4.36	3.18	3.66	3.84	4.18	3.60

With GWR analysis, the interpolated local R² surfaces for the 10 years were also averaged and displayed in Figure 11. Apparent spatial nonstationarity variations can be found in these surfaces. For the Peak NDVI (Figure 11(a)), Cumulative NDVI (Figure 11(f)), and Peak date (Figure 11(b)), the overall spatial patterns of the local R² values were similar, with a contiguous area extending from the north to central southeast exhibiting the greatest explained variance (Figure 11(a) and 11(f)). As shown in Figure 2, this area is mainly covered by grasses in mixed and tall grass prairies with relatively similar biophysical properties and mild climate conditions (Wang et al. 2013). Conversely, the lowest local R² values were concentrated in the southwestern Great Plains, which mainly consist of warm-season short grass steppe in a typical semi-arid climate. The local R² of the Peak date showed similar trends as those of the Peak NDVI, although the high values extended to the southeast end of the Plains. Local R² surfaces of the SOS, EOS, and LOS displayed an apparent north-south trend with the highest values in the north, indicating that cold-season prairies in the northern states were more sensitive to local climate variations. These spatial patterns of climate dependencies for vegetation phenology could improve our understanding to how vegetation responds to climate change in this important grassland eco-region. Finally, parts of southwestern Kansas and central Nebraska, which are underlain by the Ogallala Aquifer are subject to irrigation and this anthropogenic effect modifies phenology.

Figure 11. Ten-year mean local R² derived from the GWR models and spatial interpolation. (a) Peak NDVI as dependent variable; (b) Peak date as dependent variable; (c) SOS as dependent variable; (d) EOS as dependent variable; (e) LOS as dependent variable; and (f) Cumulative NDVI as dependent variable.

5. Conclusions

In this study, we calculated the phenological metrics for US Great Plains from 2000 to 2009, to analyze the spatial patterns of LSP. We also investigated the patterns of phenological metrics relative to monthly climatic factors. We found the Peak NDVI and Cumulative NDVI had very similar spatial patterns, with a gradient extending from a minimum in the west to a maximum in the east when crossing the Plains, a finding that parallels precipitation patterns in this region. The growing season metrics (SOS, EOS, and LOS) were determined by both the temperature and precipitation conditions, with the earliest SOS and vegetation Peak date as well as the longest LOS found in the southeast part of the Great Plains. Interestingly, the latest SOS was associated with the southwestern Plains in New Mexico where late season rainfall associated with the Arizona Monsoon proved critical. The EOS variable exhibited a more consistent pattern across the whole region because of its overwhelming dependence on temperature even though the statistical model indicated November precipitation was crucial for an extended growing season (cloud cover limiting nocturnal temperature minima).

Our regression analyses indicated early spring and fall dew points were the most determining factors for the phenological patterns. The GWR models revealed local spatial patterns of the relationships between the phenological indices and climatic variables, thereby confirming the GWR tool can be applied to model phenological indices effectively based on climatic factors. The extracted remote sensing-based phenological metrics can be used effectively in analyses of phenology-assisted crop mapping and regional biomass production assessment in this region in future studies where field-based data covering large areas are either lacking or difficult to obtain.