Land surface phenology along urban to rural gradients in the U.S. Great Plains

Remote Sensing of Environment 165 (2015) 42–52

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Land surface phenology along urban to rural gradients in the U.S. Great Plains J.J. Walker a,⁎, K.M. de Beurs b, G.M. Henebry a a b

Geospatial Sciences Center of Excellence, South Dakota State University, Brookings, SD, USA Department of Geography and Environmental Sustainability, The University of Oklahoma, Norman, OK, USA

a r t i c l e

i n f o

Article history: Received 13 June 2014 Received in revised form 29 March 2015 Accepted 15 April 2015 Available online xxxx Keywords: Urban heat island (UHI) Phenology MODIS Land surface temperature (LSP) NDVI Great Plains

a b s t r a c t The elevated surface and air temperatures of urban environments can influence the timing of vegetation growth dynamics within and across city boundaries. We examined patterns of land surface phenology (LSP) throughout the U.S. Great Plains region, which contains diverse metropolitan areas embedded within a predominately agricultural landscape. We assembled a time series (2002–2012) of Moderate Resolution Imaging Spectroradiometer (MODIS) Nadir BRDF-Adjusted Reflectance (NBAR) data and land surface temperature data at 500 m and 1000 m spatial resolution, respectively. We derived measures of the vegetated land surface and the thermal regime of the growing season at 8-day intervals using the Normalized Difference Vegetation Index (NDVI) and Accumulated Growing Degree-Days (AGDD). Fitting the convex quadratic LSP model of NDVI as a function of AGDD yielded two phenometrics – Peak NDVI and Thermal Time to Peak – and one model fit metric – the adjusted coefficient of determination (r2adj) – for each pixel per growing season. We linked the phenometrics with impervious surface area (ISA) data extracted from the U.S. Geological Survey National Land Cover Dataset (NLCD) to characterize differences in timing and amplitude of peak greenness between urban areas and their surrounding landscapes. Our results reveal the broad control of climatic conditions and moisture availability on phenological patterns across urban to rural gradients, with drier, southern cities displaying more varied responses of peak greenness timing and amplitude to urban intensity. Published by Elsevier Inc.

1. Introduction Human populations are increasingly found in urban settings. In the United States, nearly half of non-metropolitan counties lost population between 1988 and 2008 (McGranahan, Cromartie, & Wojan, 2010). Rural flight is particularly pronounced across the agricultural U.S. interior, where in-migration to metropolitan areas has increased steadily since the 1950s. At the same time, the population of rural counties has declined (Brown, Johnson, Loveland, & Theobald, 2005; Gutmann, Parton, Cunfer, & Burke, 2005; Wilson, 2009) and cropland area has expanded (Wright & Wimberly, 2013). The demographic shift from scattered rural communities to denselypopulated metropolitan areas has transformative social (Alberti, 2005), biogeochemical (Grimm et al., 2008), and ecological (McDonnell & Pickett, 1990) implications for affected populations and the surrounding environment. The capacity of urban areas to modify local climate through the urban heat island (UHI) effect has long been recognized (e.g., Howard, 1833). As an unintentional byproduct of urbanization, the UHI develops ⁎ Corresponding author at: Western Geographic Science Center, United States Geological Survey, 520 N. Park Ave., Tucson, AZ 85719, USA. Tel.: +1 520 670 5544; fax: +1 520 670 3300. E-mail address: [email protected] (J.J. Walker).

http://dx.doi.org/10.1016/j.rse.2015.04.019 0034-4257/Published by Elsevier Inc.

from the replacement of natural ground cover with impervious surface materials such as concrete, asphalt, and metal (Landsberg, 1981). Higher fractions of impervious surface area increase the amount of incoming solar radiation partitioned to sensible heat flux, in contrast to the higher proportions of latent heat flux over vegetated surfaces (Avissar, 1996). Urban geometry additionally hinders cooling by blocking sky-view and promoting the absorption of longwave radiation by adjacent surfaces (Oke, Johnson, Steyn, & Watson, 1991). The resulting radiative cooling differences produce a thermal contrast between urban and rural surroundings that can alter regional climate (Oke, 1987) and influence precipitation patterns (Dixon & Mote, 2003; Niyogi et al., 2011; Schmid & Niyogi, 2013). The existence of climatic differences between urban and rural areas is almost universal (Luo, Sun, Ge, Xu, & Zheng, 2007); the UHI has been documented in cities with fewer than 10,000 residents (Karl, Diaz, & Kukla, 1988) and in diverse locations around the world (Clinton & Gong, 2013; Peng et al., 2012). The particular thermal manifestation over each city, however, can vary significantly according to site-specific ecological and climatic drivers (Arnfield, 2003). UHI magnitude and sign are strongly influenced by meteorological conditions (Oke, 1981), the type and amount of vegetation replaced by urbanization, and the intensity and geographical extent of development (Clinton & Gong, 2013; Imhoff, Zhang, Wolfe, & Bounoua, 2010; Zhang, Imhoff, Wolfe, & Bounoua, 2010). In temperate

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

locations, air temperatures in urban areas tend to be warmer than their rural counterparts primarily at night, and have higher minimum temperatures and a constrained diurnal range (Oke, 1982). In desert environments, the low vegetative fraction in outlying areas contrasts with the irrigated landscape maintained by urban residents, setting up an inversion of the UHI that yields a heat sink during summer months (Buyantuyev & Wu, 2012). Since temperature is one of the environmental cues for springtime growth initiation in many plant species (Rathcke & Lacey, 1985; White, Thomton, & Running, 1997), the influence of urbanization on vegetation phenology may provide a preview of the consequences of future climate warming (White, Nemani, Thornton, & Running, 2002; Zhang, Friedl, Schaaf, & Strahler, 2004b; Ziska et al., 2003). Plant development theory suggests that higher temperatures should prompt earlier growth in warmth-sensitive species, as documented by the empirical observations of multiple studies (e.g., Menzel, 2000; Parmesan & Yohe, 2003; Schwartz, Ahas, & Aasa, 2006). Ground-based studies of the UHI effect on the timing of flowering and growth of individual species have both confirmed the general validity of this premise and documented its complex relationship to urbanization (e.g., Gazal et al., 2008; Luo et al., 2007; Mimet et al., 2009). Urban development in the United States is heterogeneous and patchy; green spaces and parks can be interspersed within densely settled core areas, while large industrial complexes are often found in the urban periphery. The mosaic of surfaces in urban settings (Arnfield, 2003; Cadenasso, Pickett, & Schwarz, 2007) and the variability of species-specific sensitivity to environmental growth cues (Neil, Landrum, & Wu, 2010) preclude a synchronized and predictable response of phenological timing to urbanization intensity. Mimet et al. (2009) found that the initiation of deciduous budburst stages did not linearly advance with distance along the rural–urban gradient due to the buffering influence of local environments. In a multicity study of the budburst timing of native deciduous trees, Gazal et al. (2008) concluded that although each city exhibited urban/rural phenological differences, the temperate model of earlier budburst in urban areas did not apply to tropical cities, and vegetation responses often contradicted classic plant growth paradigms. Luo et al. (2007) found that three different species exhibited variable phenological responses to air temperatures in Beijing. Remote sensing studies of the land surface phenology (LSP) (de Beurs & Henebry, 2004; Ganguly, Friedl, Tan, Zhang, & Verma, 2010; Henebry & de Beurs, 2013; Morisette et al., 2009) in and around urban areas provide synoptic views of the aggregate landscape response to climatic conditions and fractional impervious surface area. The broader spatial scale dampens the heterogeneity and variability of plant dynamics, yielding a more generalized view of composite phenological trends. As seen via remote sensing data at moderate spatial resolutions, urbanization in the forested northeastern United States accelerates the timing of vegetation development within and around the urban extent, resulting in earlier growth and longer growing seasons than in the rural perimeter (Fisher, Mustard, & Vadeboncoeur, 2006; White et al., 2002; Zhang, Friedl, Schaaf, Strahler, & Schneider, 2004a). Studies of other forested biomes have documented a similar increase in growing season length (Zhang et al., 2004b). Researchers have investigated the impact of UHIs on land surface phenology within representative bioclimatic regions (White et al., 2002) and individual cities around the world (e.g., Luo et al., 2007), but to date none have examined the influence of UHIs on the land surface phenology of the U.S. agricultural heartland. The changing nature of population distributions across the U.S. Great Plains presents an opportunity to study the emergent effect of UHIs on local and regional phenological patterns. Great Plains cities share broad commonalities of ecological setting: the rural vegetative composition is predominantly an herbaceous mixture of agriculture crops and grasslands, and the cities are sheltered from the meteorological influences of large water bodies and mountain ranges. The diversity of metropolitan sizes and densities, the range of agricultural crop species, and the broad

43

latitudinal and longitudinal extent, meanwhile, yield biophysical and climatic variations that drive multiple potential agents of phenological diversity. The U.S. Great Plains thus provide an ideal backdrop against which to investigate the spectrum of vegetation signal responses that accompany large-scale urban displacement of herbaceous land cover. In this study we explore the relationship between accumulated thermal time and indices of landscape greenness to extract measures of phenological transitions (phenophases) across the urban and peri-urban areas of the Great Plains over the past decade (2002–2012). Our specific objective is to link the phenology metrics with city-specific measures of urbanization intensity in order to characterize the influence of climate and location drivers on patterns of peak greenness.

2. Study site 2.1. U.S. Great Plains The study domain is the expansive and variously-defined (Rossum & Lavin, 2000) region of the U.S. Great Plains (Fig. 1). The climate ranges from the continental seasons of northern North Dakota, with warm, humid summers and long, cold, snowy winters, to the transitional semi-arid climate of Dallas–Fort Worth, Texas. Precipitation rates generally decrease along an east to west gradient, from a mean 20th century high of 1180 mm in the east Texas municipalities of Longview, Texarkana, and Tyler, to a low of 315 mm approximately 725 km west, in the town of Odessa, Texas (NOAA NCDC, 2014a). During the time frame of this study, the region experienced area-specific drought conditions in multiple years: 2002, 2006, 2007, 2011, and 2012 (Dong et al., 2011; NOAA NCDC, 2014b). The unusually warm spring of 2012 prompted early flowering and budding across much of the region (Ault et al., 2013). The ecological setting is broadly herbaceous, with specific agricultural crop or plant communities dictated by local climate and biophysical conditions. Warm-season crops such as maize, alfalfa, spring wheat, and soybeans predominate in the central and northern regions of the study area, while the southern latitudes are dominated by cool-season crops such as cotton, winter wheat, and sorghum. Grasslands are common along the drier, western reaches of the study area. The types and growth patterns of vegetation within metropolitan limits in the Great Plains can contrast markedly with those of the surrounding rural land cover due to irrigation and the cultivation of non-native species.

2.2. Urban areas We identified 51 cities of interest (Table 1) by extracting all urban areas from the 2010 U.S. Census Bureau data (U.S. Census Bureau, 2010) that fell within the study area and fit the general ecological setting study requirements. The U.S. Census Bureau defines urban areas as those which consist of densely developed territory containing 50,000 or more people. The population of the selected metropolitan areas ranges from 50,440 (Grand Island, Nebraska) to 5.1 million (Dallas–Fort Worth, Texas), while the geographical extent of each urban footprint varies from 5365 ha (Manhattan, Kansas) to 460,794 ha (Dallas–Fort Worth, Texas) (U.S. Census Bureau, 2010).

3. Data We derived temperature and vegetation index data from Moderate Resolution Imaging Spectroradiometer (MODIS) products. All products were acquired through the NASA internet data portal (http://reverb. echo.nasa.gov) for the five ISIN tiles that encompass the study area: h09v05, h10v05, h11v05, h10v04, and h11v04.

44

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

Fig. 1. (left) The 51 study cities in the U.S. Great Plains states of Iowa (IA); Kansas (KS); Minnesota (MN); Missouri (MO); Nebraska (NE); North Dakota (ND); Oklahoma (OK); South Dakota (SD); and Texas (TX), shown against the major 2006 National Land Cover Dataset (NLCD) land cover types. (right) The study cities shown against the mean annual precipitation of the National Climate Dataset Center (NDCD) climate divisions.

We then calculated accumulated growing degree-days (AGDDs) in each year from the GDD values:

3.1. Land surface temperature data Researchers have traditionally relied on discrete in situ measurements of air surface temperatures to study UHI phenomena. The advent of satellite data has made it possible to retrieve areal measurements of the emitted radiance from the terrestrial surface, which are translated via radiative transfer equations to land surface brightness temperatures (Weng, Lu, & Schubring, 2004). The advantages of synoptic coverage and temporal consistency have driven the use of remote sensingderived measures of land surface temperature to study UHI phenomena, given the inherent correlation between surface temperatures and surface conditions and processes (Imhoff et al., 2010; Nichol & Wong, 2005; Owen, Carlson, & Gillies, 1998; Zhang, Imhoff, Bounoua, & Wolfe, 2012) and the ability to monitor concurrent changes in the underlying causes of UHIs (higher ISA and lower vegetation fraction) (Gazal et al., 2008). Air and surface temperatures are not perfectly correlated (Arnfield, 2003; Roth, Oke, & Emery, 1989); surface temperatures are higher and more controlled by surface characteristics than coincident air temperatures (Arnfield, 2003; Nichol, 1996; Streutker, 2002). The wide variety of urban surface materials arranged in different morphological configurations yields increased variability of surface temperatures compared to air temperatures. Surface temperature UHIs (SUHIs) typically register the greatest temperature discrepancies between urban and rural areas during the day and the smallest differences at night, which is the inverse of the typical pattern of air temperature UHIs (Roth et al., 1989). We downloaded the MODIS land surface temperature (LST) product at 1 km spatial resolution and 8-day composite temporal resolution (product MOD11A2, version 5) from the first available dataset in 2002 to the last available dataset in 2012. MOD11A2 data are produced for each 8-day period by averaging two to eight days of the daily daytime and nighttime LST data, respectively (Wan, 2009). We screened cloudy or missing data from each pixel prior to calculating the simple average of day and night land surface temperatures for each 8-day time period t, converting from Kelvin to degrees Celsius, and determining the growing degree-days using a base of 0 °C: GDDt ¼ max Tday þ Tnight





=2–273:15; 0 :

ð1Þ

AGDDt ¼ AGDDt‐1 þ 8  GDDt

ð2Þ

where AGDDt is the number of growing degree-days accumulated from the beginning of the year until time period t. Multiplying by a factor of 8 accounts for the 8-day composite period of MOD11A2 data. We used a base of 0 °C since this threshold is often applied in the modeling of annual crops and perennial grasslands (de Beurs & Henebry, 2008; Henebry, 2013). We subsequently resampled the data to 500 m using cubic convolution to match the higher spatial resolution of the surface reflectance data. 3.2. Land surface reflectance data We downloaded the Collection 5 MODIS Nadir BRDF (Bidirectional Reflectance Distribution Function)-Adjusted Reflectance (NBAR) data product (MCD43A4) at 500 m resolution and 16-day temporal resolution from January 2002 to December 2012. The MCD43A4 product uses multi-angle surface reflectance values to model the data that would have been obtained from a nadir view and the mean solar zenith angle of the 16-day compositing period (Lucht, 1998; Lucht, Schaaf, & Strahler, 2000; Schaaf et al., 2002). The NBAR dataset is produced every 8 days using combined data from MODIS instruments on the Aqua and Terra satellites. For each composite image we calculated the Normalized Difference Vegetation Index (NDVI) using the NBAR red and near-infrared bands, as follows: NDVI ¼ ðNIR–RedÞ=ðNIR þ RedÞ:

ð3Þ

The individual NDVI scenes were assembled into time series of 46 images for each year. 3.3. Land cover and impervious surface data The 2006 National Land Cover Database (NLCD) (Fry et al., 2011) supplied the percent of impervious surface area (ISA) at 30 m resolution. ISA data are derived from Landsat 7 ETM + and IKONOS data

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

45

Table 1 Characteristics of the 51 study area cities: population and land extent encompassed by the urban area boundary (U.S. Census Bureau, 2010); latitude, longitude, and 20th century mean annual precipitation (NOAA NCDC, 2014a) of the encompassing climate division; and percent of the pixels within the urban boundary for which the Land Surface Phenology (LSP) model consistently returned r2 values greater than 0.8.

Abilene, TX Amarillo, TX Ames, IA Bismarck, ND Cape Girardeau, MO Cedar Rapids, IA Columbia, MO Dallas–Fort Worth, TX Davenport, IA Denton, TX Des Moines, IA Dubuque, IA Fargo, ND Grand Forks, ND Grand Island, NE Iowa City, IA Jefferson City, MO Joplin, MO Kansas City, MO Lawrence, KS Lawton, OK Lee's Summit, MO Lincoln, NE Longview, TX Lubbock, TX Manhattan, KS Mankato, MN McKinney, TX Midland, TX Minneapolis−St. Paul, MN Norman, OK Odessa, TX Oklahoma City, OK Omaha, NE Rochester, MN San Angelo, TX Sherman, TX Sioux City, IA Sioux Falls, SD Springfield, MO St. Cloud, MN St. Joseph, MO St. Louis, MO Texarkana, TX Topeka, KS Tulsa, OK Tyler, TX Waco, TX Waterloo, IA Wichita Falls, TX Wichita, KS

City population

Land area (km2)

Latitude (°N)

Longitude (°W)

Precipitation (mm)

LSP model (%)

110,421 196,651 60,438 81,955 52,900 177,844 124,748 5,121,892 280,051 366,174 450,070 67,818 176,676 61,270 50,440 106,621 58,533 82,775 1,519,417 88,053 94,457 85,081 258,719 98,884 237,356 54,622 57,584 170,030 117,807 2,650,890 103,898 126,405 861,505 725,008 107,677 92,984 61,900 106,494 156,777 273,724 110,621 81,176 2,150,706 78,162 150,003 655,479 130,247 172,378 113,418 99,437 472,870

141.8 210.3 60.2 100.4 90.7 216.4 160.1 4607.9 358.0 376.0 519.5 87.5 182.0 63.3 72.9 118.0 103.6 166.8 1755.6 78.8 114.1 107.0 229.1 215.1 249.8 53.6 68.3 191.9 136.9 2646.5 116.2 152.6 1063.5 702.4 131.0 120.9 93.0 140.8 166.2 368.6 130.1 108.4 2392.2 166.8 206.1 870.0 233.7 233.7 161.2 130.4 556.2

32.428 35.198 42.032 46.820 37.333 41.999 38.956 32.813 41.517 33.093 41.619 42.506 46.862 47.906 40.926 41.681 38.576 37.101 39.040 38.956 34.618 38.898 40.807 32.517 33.559 39.194 44.167 33.205 31.997 44.978 35.214 31.865 35.505 41.235 44.020 31.441 33.683 42.482 43.533 37.149 45.574 39.757 38.631 33.439 39.042 36.105 32.310 31.542 42.496 33.903 37.680

99.747 101.849 93.634 100.800 89.581 91.663 92.328 96.972 90.513 97.066 93.656 90.691 96.821 97.057 98.355 91.559 92.203 94.499 94.595 95.259 98.416 94.376 96.680 94.779 101.893 96.593 93.991 96.677 102.107 93.281 97.445 102.401 97.524 96.035 92.480 100.460 96.583 96.402 96.736 93.287 94.195 94.830 90.341 94.071 95.692 95.901 95.302 97.166 92.384 98.519 97.327

597 479 827 416 1106 861 960 869 861 869 827 843 557 471 601 861 1035 1095 938 917 711 938 723 1180 479 877 736 869 479 717 874 315 874 723 781 641 869 764 580 1095 675 938 960 1180 917 1016 1180 869 843 597 687

93.4 99.0 100 99.6 100 99.6 99.9 99.1 99.3 99.6 99.9 100 99.7 99.4 100 99.8 99.8 99.9 99.8 100 100 100 99.8 99.8 94.6 100 98.2 99.1 63.4 99.0 100 46.0 99.1 99.9 100 80.3 100 99.9 100 99.9 99.8 100 99.8 100 99.6 99.9 100 99.8 100 98.8 99.8

and discriminate manmade surfaces (i.e., asphalt, concrete) from natural ones (Yang, Huang, Homer, Wylie, & Coan, 2003). We computed the percentage of ISA contained within each 500 m MODIS pixel in the study area (Fig. 2). The continental-scale datasets were subset to the dimensions of the study domain to minimize processing load.

can consistently and reliably characterize herbaceous land cover dynamics in temperate regions (Brown & de Beurs, 2008; de Beurs & Henebry, 2004, 2005, 2008, 2010). The quadratic regression model requires the estimation of just three parameters, each of which has a straightforward ecological interpretation (de Beurs & Henebry, 2005). The convex quadratic (CxQ) model takes the form:

4. Methods

NDVI ¼ α þ β  AGDD ‐γ  AGDD2

4.1. Modeling of AGDD as a function of NDVI

where AGDD is the accumulated growing degree-days (°C); α is the background NDVI value at the start of the growing season; and the parameters β and γ describe the slope and curvature of the annual green-up and senescence period. The peak of modeled NDVI curve can be expressed in terms of the fitted parameter coefficients:

Researchers have investigated diverse approaches and techniques for extracting phenological benchmarks from seasonal vegetation growth curves (e.g., de Beurs & Henebry, 2010). Particular interest has centered on fitting modeled relationships of NDVI dynamics derived from satellite imagery (Jönsson & Eklundh, 2002; Zhang et al., 2004b). Modeling the relationship of NDVI to AGDD as a quadratic function

Peak NDVI ¼ α‐β2  ð4  γÞ‐1 :

ð4Þ

ð5Þ

46

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

deleted pixels from the study if their associated metrics were beyond reasonable phenological limits, e.g., pixels with peak NDVI dates prior to DOY 30 or after DOY 330. On average, 1.4% of each city's pixels were removed. In each city we calculated the mean timing and height of peak NDVI within impervious surface area (ISA) intervals 0–15%, 16–25%, 26–45%, 46– 65%, and 66–100%, which represent a gradient of rural to most intensively urbanized pixels across the landscape. 4.3. Statistical analysis All descriptive statistics and regression analyses were performed using the R statistical package (v. 2.15.2). Image processing and vector-based analyses were performed using ENVI (v. 5.0) and ESRI's ArcGIS (v. 10.2). 5. Results 5.1. Land surface phenology model

To avoid the retrieval of unreasonable seasonal profiles, we set model tolerances for the season length (30 to 300 days); the alpha parameter (−100 to 0.70); peak position (20 to 6000); peak height (0.1 to 1); and adjusted r2 (0.8 to 1). A fitted model was accepted only if all parameter requirements were satisfied. We fitted the CxQ model to each pixel time series separately for every year to derive annual phenology metrics of the peak NDVI height and the thermal time to peak height in AGDD; we additionally retrieved the coefficient of determination (r2) that describes the goodness-of-fit of the model to the data. We subsequently calculated the Day of Year (DOY) equivalent of the thermal time to peak height via linear interpolation. We specifically focused on the peak of the growing season because it is readily identifiable through the fitted coefficients and is less likely to be affected by noise than the start of the growing season, resulting in a more stable estimation (de Beurs & Henebry, 2010; Henebry & de Beurs, 2013).

The severe drought that affected the southern portion of our study area in 2011 and the anomalously warm early months of 2012 resulted in higher numbers of model anomalies, particularly in the arid southwest. Since our primary interest is in the model performance during typical climatic conditions and in more temperate environments, we report model results over the entire study site for the 2002–2010 time period only. The map of average growing season peak NDVI values (Fig. 3, left) clearly shows lower amplitudes in urban areas, which is expected since urban environments in this region are associated with lower vegetation density. The map of average thermal times to peak (Fig. 3, right) shows that most urban areas exhibit peaks at higher AGDD than their surroundings; this phenomenon is particularly evident in cities in Oklahoma and eastern Texas. Since the thermal time to peak NDVI depends on the individual rate of heat accumulation in each pixel, the translation of a delayed thermal peak to Day of Year does not necessarily result in the same relative relationship of peak timing between an urban area and its surroundings. The maps of average adjusted coefficients of determination (r2adj) and the number of years the CxQ model failed the selection criteria (Fig. 4) demonstrate the robustness of the model under most conditions and in all but the most arid regions of the study region. All cities show wellfitting models in each year for more than 95% of the pixels, except for the dry west Texas cities of Abilene, Lubbock, Midland, Odessa, and San Angelo (Table 1).

4.2. Delineation of urban extent

5.2. Peak NDVI timing vs. mean %ISA

The lack of a universally-accepted delineation of “urban space” complicates attempts to track signals along gradations of urban development (Montgomery, Stren, Cohen, & Reed, 2003; Potere, Schneider, Angel, & Civco, 2009; Raciti, Hutyra, Rao, & Finzi, 2012). The mixture of natural and man-made surfaces in urban areas undermines the validity of associating gradations in phenological metrics with a simple measure of linear distance from an arbitrary urban center (e.g., Mimet et al., 2009). A more consistent method of tracking the influence of urban development on phenology incorporates the proportion of man-made surface as a proxy for urbanization levels, since the amount of impervious surface area correlates with the magnitude of the urban heat island effect (Imhoff et al., 2010; Roth et al., 1989; Zhang et al., 2012). We restricted our rural analysis to an area within 40 km of the Census urban boundary of each city to fully encompass urban, suburban, and rural extents. Each 500 m pixel within this defined geographical area was linked with phenology and urbanization intensity metrics. We removed pixels 10 km beyond the metropolitan boundary with a mean ISA greater than 5% to avoid biasing the rural signal with pixels associated with nearby towns or isolated, highly developed areas. We additionally

The 51 study cities display individualized relationships between the DOY of peak NDVI and the %ISA in each MODIS pixel (Fig. 5). Nonetheless, the cities can be loosely grouped into categories that share similar profile characteristics dictated by climate and land cover commonalities. Cities in and around the western Corn Belt (here Iowa, Nebraska, Kansas, and Missouri) typically reveal a negative relation between %ISA and average peak timing (DOY) for 2002–2010, resulting in earlier growing season peaks in pixels with higher %ISA. Cities in the southern portions of the study region typically have a positive relation between %ISA and peak timing (DOY); i.e., higher %ISA pixels are associated with later growing season peaks. Cities far north and in the transition zone between these two regions reveal no significant relationships between %ISA and the timing of growing season peak. The unprecedented false spring of 2012 (Ault et al., 2013) prompted the earliest peak dates across the urban/rural gradient in many of the eastern cities; in contrast, the drier cities in the south register earliest peaks in either 2003 or 2011, while the earliest peak in the only dry northern city (Bismarck, North Dakota) was in 2006. In general, the early onset of warmth led to a more urban-mediated response, with a wider range of

Fig. 2. Percent of impervious surface area (ISA) in and around the metropolitan area of Omaha, NE-Council Bluffs, IA. The external buffer extending to 40 km beyond the census boundary is not depicted here.

The quantity of AGDD or “thermal time” required to reach the peak NDVI in each pixel can also be computed from the fitted parameter coefficients: Thermal Time to Peak ¼ −β  ð2  γÞ−1 :

ð6Þ

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

47

Fig. 3. Average (2002–2010) land surface phenology model output showing (left) height of the peak NDVI, and (right) thermal time to peak NDVI. Urban areas appear as outlined polygons. Contrasts between the urban areas and their surrounding landscapes are evident for most cities.

days between the earliest and latest dates of peak NDVI across ISA categories (Fig. 6). 5.3. Variability of peak NDVI timing vs. mean %ISA along latitudinal and longitudinal extents We examined the variability (i.e., the standard deviation [SD]) in peak NDVI timing as a function of longitude and latitude and %ISA to determine whether urbanization intensity was associated with constrained peak

greenness measures across the region. For ease of analysis and illustration purposes, we grouped the cities by longitude to align roughly with high (east of 96° W), medium (between 96° and 98° W), and low (west of 98° W) regional precipitation patterns. The resulting distributions clearly show the controlling effect of latitude and longitude on urban phenological patterns (Fig. 7). At every ISA level, annual phenological variability tends to increase with more southerly latitudes and more westerly longitudes. We used analysis of covariance (ANCOVA) to analyze whether the increase in variability with decreasing latitude changes across the east–

Fig. 4. (left) Average adjusted coefficients of determination (r2adj) for the land surface phenology models (2002–2010) and (right) the number of years that the model failed for each pixel in the study region. Overall, the quadratic regression models displayed good fits with most pixels showing an average r2adj greater than 0.8. In the arid southwestern part of the study region the thermal time model frequently fails to fit.

48

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

Fig. 5. (left) Pearson correlations (r) between the mean timing of peak NDVI (2002–2010) and percent of impervious surface area (%ISA) categories: 0–15%; 16–25%; 26–45%; 46–65%; and 66–100%. Data were extracted from all pixels within U.S. Census urban boundaries and surrounding 40 km extents. Relationships were assessed as strongly (p b 0.01), moderately (0.01 ≤ p b 0.05), or weakly (0.05 ≤ p b 0.10) significant; not significant (0.10 ≤ p b 0.20); and no relation (p ≥ 0.20). (right) The metropolitan areas of Oklahoma City, OK, and Dubuque, IA, display relationships characteristic of most southern cities (later peak timing associated with higher % ISA) and mid-latitude cities (earlier peak timing associated with higher %ISA). The background image shows average annual temperatures (PRISM, 2015) indicating relatively warmer areas in the south and cooler areas in the north.

west extent of our study area for each of the ISA percentages of interest. ANCOVA tests the null hypothesis of no difference between the slopes of the regression lines fitted to the data of each longitude category. Results suggest that the central and eastern longitudinal categories are significantly different from the western category for ISA ranges of 25% and 50% (p b 0.05). No such interaction exists for ISA 75%, meaning that the magnitude and rate of change of the variability in the three categories were statistically equivalent from east to west along the latitudinal extent.

5.4. Peak NDVI height vs. mean %ISA The height of peak NDVI displays far less variability in individual city responses (Fig. 8). All cities display a negative relationship of peak NDVI

height to %ISA that is relatively immune to climatic conditions across time. Deviations are seen only in the semi-arid, southwestern cities, where the discrepancies between the categories uniformly decrease in drought years (i.e., 2006, 2011) to the levels of the most urbanized pixels.

6. Discussion The Great Plains region encompasses a broad latitudinal extent and a pronounced east–west precipitation gradient that exhibits considerable year to year variation in timing and amount. Our results highlight the overarching influence of these factors—i.e., growing season length and aridity—on the phenological behavior exhibited by the metropolitan

Fig. 6. Boxplots of the number of days between the earliest and latest dates of annual peak NDVI for all 51 study cities. For each city the average peak timing was calculated within five impervious surface area (ISA) intervals: 0–15%; 16–25%; 26–45%; 46–65%; and 66–100%. The solid horizontal bar across each box indicates the median value of the range; the lower and upper edges of the box represent the first and third quartiles, respectively; and the end of each whisker shows the highest or lowest value within 1.5 times the inter-quartile range.

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

49

Fig. 7. Standard deviation of the day of year (DOY) of peak NDVI in pixels containing (top) 0–15%, (middle) 25–45%, or (bottom) 65–100% impervious surface area (ISA) in all study area cities from 2002 to 2012. All pixels within each urban area and its surrounding 40 km extent were used in the analysis.

areas of this study. In general, the relative location of each city along the latitudinal and longitudinal axes of the Great Plains study domain is a more reliable predictor of broad phenology patterns across the rural to urban gradient than characteristics such as development intensity. Climatic conditions and moisture availability dictate the agricultural crops and natural vegetation species found in the peri-urban and rural

environs; the phenologies of these land cover types are likely to be relatively uniform in areas subject to similar climatic and environmental conditions. The variety of peak NDVI timing and height patterns from cities within common climate divisions shown in Fig. 1, meanwhile, demonstrates the complex influence of secondary factors in controlling site-

50

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

Fig. 9. Peak NDVI timing plotted as a function of year (2002–2012) within impervious surface area (ISA) intervals 0–15% and 66–100% for three Texas cities that span the longitudinal range of annual precipitation: Lubbock (479 mm), Sherman (869 mm), and Texarkana (1180 mm). Error bars indicate one standard deviation.

Fig. 8. Peak NDVI plotted as a function of year (2002–2012) within impervious surface area (ISA) intervals 0–15% and 66–100%. The profiles display the responses typical of high- and mid-latitude cities (Rochester, MN, and Ames, IA); lower-latitude eastern (wetter) cities (Dallas–Ft. Worth, TX); and (d) lower-latitude western (semi-arid) cities (San Angelo, TX). Error bars indicate one standard deviation.

specific responses to mean and anomalous weather regimes on a seasonal basis. These nested patterns allow us to summarize regional trends on a north/south and east/west basis. Latitude controls the length of day and affects the amount of insolation received during the growing season, which in turn determines the potential variation in peak NDVI timing. Northern Great Plains cities show a relatively more consistent and constrained range of peak NDVI dictated by the relatively shorter growing season, as peak NDVI timing cannot deviate substantially from the narrow range of favorable growing conditions. Southern cities display NDVI timing profiles indicative of a longer growing season in which vegetation can achieve peak growth in early spring or late fall depending on prevailing weather conditions. The longitudinally-influenced precipitation amount is another critical determinant of NDVI peak timing and peak height patterns. Precipitation has been identified as the most important factor controlling Great Plains grasslands production (Sala, Parton, Joyce, & Lauenroth, 1988; Tieszen, Reed, Bliss, Wylie, & DeJong, 1997). Drier Texas cities display more annual and inter-annual variation and less synchronized patterns of timing among the different ISA categories (Fig. 9). Abnormal temperature and drought conditions affect the patterns of peak growth differently across the region. The drought experienced by much of the region in 2006 (Dong et al., 2011) generally resulted in a contraction of the ISA-modulated amplitude response in the western Texas cities, where the rural vegetation signal was relatively indistinguishable from the urban one during periods of drought (Fig. 8). The premise of most UHI-related urban phenology studies is that observed differences in the timing of vegetative growth events are attributable to the effect of the urban environment, which alters the timing of growth initiation and peak. However, land cover differences between the urban environment and the surroundings influence the timing of

growth due to the different species-specific requirements. Within arid cities, urban vegetation can become unsynchronized with the phenological signals of the surrounding natural desert vegetation due to the widespread use of irrigation within the city and the presence of nonnative species (Buyantuyev & Wu, 2012). Since this study did not control for differences in land cover composition within and surrounding the urban area, we cannot definitively state that the phenological differences are due to the UHI-prompted earlier start to the growing season. However, given the range of vegetation communities covered by the NLCD land cover classes and the rapid pace of land cover conversion in urbanizing areas, it is doubtful that controlling for land cover composition would be effective at isolating the signal. 7. Conclusion This study documents the complex relationship between urban intensity and phenology across a region with broadly similar environmental characteristics. Urban landscape structure and form, speciesspecific vegetative differences, ecological context, and seasonal climate variations shape highly individualized profiles of each city's phenological response across the urban to rural gradient. On a local level, the variety of responses complicates the critical incorporation of vegetation effects into models of urban surface boundary layer dynamics at relevant scales (102–104 m) (Grimmond et al., 2010, 2011). On a more general level, the variation underscores the inadvisability of extrapolating observations about the influence of urbanization on vegetative growth from individual case studies to regional extents. The patterns that emerged in this study demonstrate the general control of climatic and geographic context on phenological response across urban gradients, and suggest that the inclusion of finer-scale data and additional urban parameters is necessary for attributing and interpreting the distinctive features of each city's phenological signature. Future research will consider physical characteristics of the urban landscape that were not explored in this analysis—e.g., the arrangement of green spaces, fragmentation of the urban fringe, and compactness of the urban core—for their ability to contribute more nuanced insights into the underlying causes of the effect of built-up environments on vegetation response. An additional avenue for exploration is the spatial and temporal resolution of the remote sensing data used for the extraction of vegetated and urban surfaces. The high temporal resolution of

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

the MODIS data used in the study is offset by its relatively coarse spatial resolution, which subsumes considerable detail of the heterogeneous urban areas. The spatial resolution was appropriate for our study objective of mapping phenology patterns across broad geographical extents. The free availability of images in the Landsat archive and the development of synthesis products such as the Web-Enabled Landsat Data (WELD) project (Roy et al., 2010) hold promise for basing future urban phenology analyses on time series that simultaneously provide details of plant development as well as surface materials and composition.

Acknowledgments Research is supported by NASA IDS project NNX12AM89G: Storms, Forms, and Complexity of the Urban Canopy: How Land Use, Settlement Patterns, and the Shape of Cities Influence Severe Weather. We would like to thank Paul de Beurs for the development of the land surface phenology derivation software used in this project. We also appreciate the insightful comments and suggestions provided by two anonymous reviewers.

References Alberti, M. (2005). The effects of urban patterns on ecosystem function. International Regional Science Review, 28(2), 168–192. http://dx.doi.org/10.1177/0160017605275160. Arnfield, A.J. (2003). Two decades of urban climate research: A review of turbulence, exchanges of energy and water, and the urban heat island. International Journal of Climatology, 26, 1–26. http://dx.doi.org/10.1002/joc.859. Ault, T.R., Henebry, G.M., de Beurs, K.M., Schwartz, M.D., Betancourt, J.L., & Moore, D. (2013). The false spring of 2012, earliest in North American record. Eos, Transactions American Geophysical Union, 94(20), 181–182. http://dx.doi.org/10.1002/2013EO200001. Avissar, R. (1996). Potential effects of vegetation on the urban thermal environment. Atmospheric Environment, 30, 437–448. Brown, M.E., & de Beurs, K.M. (2008). Evaluation of multi-sensor semi-arid crop season parameters based on NDVI and rainfall. Remote Sensing of Environment, 112, 2261–2271. http://dx.doi.org/10.1016/j.rse.2007.10.008. Brown, D.G., Johnson, K.M., Loveland, T.R., & Theobald, D.M. (2005). Rural land-use trends in the conterminous United States, 1950–2000. Ecological Applications, 15, 1851–1863. Buyantuyev, A., & Wu, J. (2012). Urbanization diversifies land surface phenology in arid environments: Interactions among vegetation, climatic variation, and land use pattern in the Phoenix metropolitan region, USA. Landscape and Urban Planning, 105(1-2), 149–159. http://dx.doi.org/10.1016/j.landurbplan.2011.12.013. Cadenasso, M.L., Pickett, S.T.A., & Schwarz, K. (2007). Spatial heterogeneity in urban ecosystems: Reconceptualizing land cover and a framework for classification. Frontiers in Ecology and the Environment, 5(2), 80–88. Clinton, N., & Gong, P. (2013). MODIS detected surface urban heat islands and sinks: Global locations and controls. Remote Sensing of Environment, 134, 294–304. http://dx.doi.org/ 10.1016/j.rse.2013.03.008. 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. http://dx.doi.org/10.1016/j.rse.2003.11. 006. de Beurs, K.M., & Henebry, G.M. (2005). A statistical framework for the analysis of long image time series. International Journal of Remote Sensing, 26, 1551–1573. http://dx. doi.org/10.1080/01431160512331326657. de Beurs, K.M., & Henebry, G.M. (2008). Northern annular mode effects on the land surface phenologies of northern Eurasia. Journal of Climate, 21, 4257–4279. http://dx. doi.org/10.1175/2008JCLI2074.1. de Beurs, K.M., & Henebry, G.M. (2010). A land surface phenology assessment of the northern polar regions using MODIS reflectance time series. Canadian Journal of Remote Sensing, 36(August), S87–S110. Dixon, P.G., & Mote, T.L. (2003). Patterns and causes of Atlanta's urban heat island– initiated precipitation. Journal of Applied Meteorology, 42, 1273–1284. Dong, X., Xi, B., Kennedy, A., Feng, Z., Entin, J.K., Houser, P.R., et al. (2011). Investigation of the 2006 drought and 2007 flood extremes at the Southern Great Plains through an integrative analysis of observations. Journal of Geophysical Research, 116(D3), D03204. http://dx.doi.org/10.1029/2010JD014776. Fisher, J.I., Mustard, J.F., & Vadeboncoeur, M.A. (2006). Green leaf phenology at Landsat resolution: Scaling from the field to the satellite. Remote Sensing of Environment, 100, 265–279. Fry, J.A., Xian, G., Homer, C.G., Yang, L., Barnes, C.A., Herold, N.D., et al. (2011). Completion of the 2006 National Land Cover Database for the conterminous United States. Photogrammetric Engineering & Remote Sensing, 77, 858–864. Ganguly, S., Friedl, M.A., Tan, B., Zhang, X.Y., & Verma, M. (2010). Land surface phenology from MODIS: Characterization of the Collection 5 global land cover dynamics product. Remote Sensing of Environment, 114, 1805–1816.

51

Gazal, R., White, M.A., Gillies, R., Rodemaker, E., Sparrow, E., & Gordon, L. (2008). GLOBE students, teachers, and scientists demonstrate variable differences between urban and rural leaf phenology. Global Change Biology, 14, 1568–1580. http://dx.doi.org/ 10.1111/j.1365-2486.2008.01602.x. Grimm, N.B., Faeth, S.H., Golubiewski, N.E., Redman, C.L., Wu, J., Bai, X., et al. (2008). Global change and the ecology of cities. Science, 319(5864), 756–760. http://dx. doi.org/10.1126/science.1150195. Grimmond, C.S.B., Blackett, M., Best, M.J., Baik, J. -J., Belcher, S.E., Beringer, J., et al. (2011). Initial results from Phase 2 of the international urban energy balance model comparison. International Journal of Climatology, 31, 244–272. http://dx. doi.org/10.1002/joc.2227. Grimmond, C.S.B., Blackett, M., Best, M.J., Barlow, J., Baik, J. -J., Belcher, S.E., et al. (2010). The international urban energy balance models comparison project : First results from phase 1. Journal of Applied Meteorology and Climatology, 49, 1268–1292. Gutmann, M.P., Parton, W.J., Cunfer, G., & Burke, I.C. (2005). Population and environment in the U.S. Great Plains. In B. Entwisle, & P.C. Stern (Eds.), Population, Land Use, and Environment: Research Directions (pp. 84–105). Washington, D.C.: The National Academies Press. Henebry, G.M. (2013). Phenologies of North American grasslands and grasses. In M.D. Schwartz (Ed.), Phenology: An Integrative Environmental Science (pp. 197–210) (2nd Ed.). New York: Springer. http://dx.doi.org/10.1007/978-94-007-6925-0-11. Henebry, G.M., & de Beurs, K.M. (2013). Remote sensing of land surface phenology: A prospectus. In M.D. Schwartz (Ed.), Phenology: An Integrative Environmental Science (pp. 385–411) (2nd Ed.). New York: Springer. http://dx.doi.org/10. 1007/978-94-007-6925-0_21. Howard, L. (1833). The Climate of London: Deduced From Meteorological Observations Made in the Metropolis and Various Places Around It. London: Harvey and Darton. Imhoff, M.L., Zhang, P., Wolfe, R.E., & Bounoua, L. (2010). Remote sensing of the urban heat island effect across biomes in the continental USA. Remote Sensing of Environment, 114, 504–513. http://dx.doi.org/10.1016/j.rse.2009.10.008. Jönsson, P., & Eklundh, L. (2002). Seasonality extraction by function fitting to time-series of satellite sensor data. IEEE Transactions on Geoscience and Remote Sensing, 40, 1824–1832. http://dx.doi.org/10.1109/TGRS.2002.802519. Karl, T.R., Diaz, H.F., & Kukla, G. (1988). Urbanization: Its detection and effect in the United States climate record. Journal of Climate, 1, 1099–1123. Landsberg, H.E. (1981). The Urban Climate. New York: Academic Press. Lucht, W. (1998). Expected retrieval accuracies of bidirectional reflectance angular sampling. Journal of Geophysical Research, 103(98), 8763–8778. Lucht, W., Schaaf, C.B., & Strahler, A.H. (2000). An algorithm for the retrieval of albedo from space using semiempirical BRDF models. IEEE Transactions on Geoscience and Remote Sensing, 38, 977–998. Luo, Z., Sun, O.J., Ge, Q., Xu, W., & Zheng, J. (2007). Phenological responses of plants to climate change in an urban environment. Ecological Research, 22, 507–514. http:// dx.doi.org/10.1007/s11284-006-0044-6. McDonnell, M.J., & Pickett, S.T.A. (1990). Ecosystem structure and function along gradients: an unexploited urban–rural opportunity for ecology. Ecology, 71, 1232–1237. McGranahan, D., Cromartie, J., & Wojan, T. (2010). Nonmetropolitan Outmigration Counties: Some are Poor, Many are Prosperous. ERR-107, U.S. Dept. of Agriculture, Economic Research Service. Menzel, A. (2000). Trends in phenological phases in Europe between 1951 and 1996. International Journal of Biometeorology, 44(2), 76–81. Mimet, A., Pellissier, V., Quénol, H., Aguejdad, R., Dubreuil, V., & Rozé, F. (2009). Urbanisation induces early flowering: Evidence from Platanus acerifolia and Prunus cerasus. International Journal of Biometeorology, 53(3), 287–298. http://dx.doi.org/10.1007/ s00484-009-0214-7. Montgomery, M., Stren, R., Cohen, B., & Reed, H. (2003). Cities Transformed: Demographic Change and Its Implications in the Developing World. Washington D.C.: National Academies Press. Morisette, J.T., Richardson, A.D., Knapp, A.K., Fisher, J.I., Graham, E.A., Abatzoglou, J., et al. (2009). Tracking the rhythm of the seasons in the face of global change: Phenological research in the 21st century. Frontiers in Ecology and the Environment, 7(5), 253–260. http://dx.doi.org/10.1890/070217. Neil, K.L., Landrum, L., & Wu, J. (2010). Effects of urbanization on flowering phenology in the metropolitan phoenix region of USA: Findings from herbarium records. Journal of Arid Environments, 74(4), 440–444. http://dx.doi.org/10.1016/ j.jaridenv.2009.10.010. Nichol, J.E. (1996). High-resolution surface temperature patterns related to urban morphology in a tropical city: A satellite-based study. Journal of Applied Meteorology, 35, 135–146. Nichol, J., & Wong, W.S. (2005). Modeling urban environmental quality in a tropical city. Landscape and Urban Planning, 73(1), 49–58. Niyogi, D., Pyle, P., Lei, M., Arya, S.P., Kishtawal, C.M., Shepherd, M., et al. (2011). Urban modification of thunderstorms: An observational storm climatology and model case study for the Indianapolis urban region. Journal of Applied Meteorology and Climatology, 50, 1129–1144. http://dx.doi.org/10.1175/2010JAMC1836.1. NOAA NCDC (2014a). National Oceanic and Atmospheric Administration (NOAA) National Climate Data Center (NCDC) Climatological Rankings. Retrieved from www.ncdc. noaa.gov/temp-and-precip/ranks.php NOAA NCDC (2014b). National Oceanic and Atmospheric Administration (NOAA) National Climate Data Center (NCDC) Climate Indices. Retrieved from www7.ncdc.noaa.gov/ CDO/CDODivisionalSelect.jsp# Oke, T.R. (1981). Canyon geometry and the nocturnal heat island: Comparison of scale model and field observation. Journal of Climatology, 1, 237–254. Oke, T.R. (1982). The energetic basis of the urban heat island. Quarterly Journal of the Royal Meteorological Society, 108(455), 1–24. Oke, T.R. (1987). Boundary Layer Climates (2nd ed.). London: Methuen.

52

J.J. Walker et al. / Remote Sensing of Environment 165 (2015) 42–52

Oke, T.R., Johnson, G.T., Steyn, D.G., & Watson, I.D. (1991). Simulation of nocturnal surface urban heat islands under ‘ideal’ conditions: Part 2. Diagnosis of causation. BoundaryLayer Meteorology, 56, 339–358. Owen, T.W., Carlson, T.N., & Gillies, R.R. (1998). An assessment of satellite remotelysensed land cover parameters in quantitatively describing the climatic effect of urbanization. International Journal of Remote Sensing, 19, 1663–1681. Parmesan, C., & Yohe, G. (2003). A globally coherent fingerprint of climate change impacts across natural systems. Nature, 421(6918), 37–42. http://dx.doi.org/10.1038/ nature01286. Peng, S., Piao, S., Ciais, P., Friedlingstein, P., Ottle, C., Bréon, F. -M., et al. (2012). Surface urban heat island across 419 global big cities. Environmental Science & Technology, 46, 696–703. http://dx.doi.org/10.1021/es2030438. Potere, D., Schneider, A., Angel, S., & Civco, D. (2009). Mapping urban areas on a global scale: Which of the eight maps now available is more accurate? International Journal of Remote Sensing, 30(24), 6531–6558. http://dx.doi.org/10.1080/01431160903121134. PRISM Climate Group (). Oregon State University. http://prism.oregonstate.edu (created 10 Jan 2015) Raciti, S.M., Hutyra, L.R., Rao, P., & Finzi, A.C. (2012). Inconsistent definitions of “urban” result in different conclusions about the size of urban carbon and nitrogen stocks. Ecological Applications, 22, 1015–1035. Rathcke, B., & Lacey, E.P. (1985). Phenological patterns of terrestrial plants. Annual Review of Ecology and Systematics, 16, 179–214. Rossum, S., & Lavin, S. (2000). Where are the Great Plains? A cartographic analysis. The Professional Geographer, 52, 543–552. http://dx.doi.org/10.1111/0033-0124.00245. Roth, M., Oke, T.R., & Emery, W.J. (1989). Satellite-derived urban heat islands from three coastal cities and the utilization of such data in urban climatology. International Journal of Remote Sensing, 10, 1699–1720. http://dx.doi.org/10.1080/01431168908904002. Roy, D.P., Ju, J., Kline, K., Scaramuzza, P.L., Kovalskyy, V., Hansen, M.C., et al. (2010). Webenabled Landsat Data (WELD): Landsat ETM+ composited mosaics of the conterminous United States. Remote Sensing of Environment, 114, 35–49. Sala, O.E., Parton, W.J., Joyce, L.A., & Lauenroth, W.K. (1988). Primary production of the central grassland region of the United States. Ecology, 69(1), 40–45. http://dx.doi. org/10.2307/1943158. Schaaf, C.B., Gao, F., Strahler, A.H., Lucht, W., Li, X., Tsang, T., et al. (2002). First operational BRDF, albedo nadir reflectance products from MODIS. Remote Sensing of Environment, 83, 135–148. Schmid, P.E., & Niyogi, D. (2013). Impact of city size on precipitation-modifying potential. Geophysical Research Letters, 40, 5263–5267. http://dx.doi.org/10.1002/grl.50656. Schwartz, M.D., Ahas, R., & Aasa, A. (2006). Onset of spring starting earlier across the Northern Hemisphere. Global Change Biology, 12(2), 343–351. http://dx.doi.org/10. 1111/j.1365-2486.2005.01097. Streutker, D.R. (2002). A remote sensing study of the urban heat island of Houston, Texas. International Journal of Remote Sensing, 23, 2595–2608. http://dx.doi.org/10.1080/ 01431160110115023.

Tieszen, L.L., Reed, B.C., Bliss, N.B., Wylie, B.K., & DeJong, D.D. (1997). NDVI, C3 and C4 production, and distributions in Great Plains grassland land cover classes. Landscape Parameterization, 7, 59–78. U.S. Census Bureau (2010). 2010 Census Urban and Rural Classification and Urban Area Criteria. Retrieved from http://www.census.gov/geo/reference/ua/urban-rural2010.html Wan, Z. (2009). Collection-5 MODIS Land Surface Temperature Products Users' Guide. Retrieved from www.icess.ucsb.edu/modis/LstUsrGuide/MODIS_LST_products_ Users_guide_C5.pdf Weng, Q., Lu, D., & Schubring, J. (2004). Estimation of land surface temperature–vegetation abundance relationship for urban heat island studies. Remote Sensing of Environment, 89(4), 467–483. http://dx.doi.org/10.1016/j.rse.2003.11.005. White, M.A., Nemani, R.R., Thornton, P.E., & Running, S.W. (2002). Satellite evidence of phenological differences between urbanized and rural areas of the eastern United States deciduous broadleaf forest. Ecosystems, 5(3), 260–273. http://dx.doi.org/10. 1007/s10021-001-0070-8. White, M.A., Thomton, P.E., & Running, S.W. (1997). A continental phenology model for monitoring vegetation resonses to interannual climatic variability. Global Biogeochemical Cycles, 11, 217–234. Wilson, S.G. (2009). Population dynamics of the Great Plains: 1950 to 2007. U.S. Census Bureau Current, Population Report P25–1137. Wright, C.K., & Wimberly, M.C. (2013). Recent land use change in the Western Corn Belt threatens grasslands and wetlands. Proceedings of the National Academy of Sciences of the United States of America, 110, 4134–4139. http://dx.doi.org/10.1073/pnas. 1215404110. Yang, L., Huang, C., Homer, C.G., Wylie, B.K., & Coan, M.J. (2003). An approach for mapping large-area impervious surfaces: Synergistic use of Landsat-7 ETM+ and high spatial resolution imagery. Canadian Journal of Remote Sensing, 29(2), 230–240. Zhang, X.Y., Friedl, M.A., Schaaf, C.B., & Strahler, A.H. (2004a). Climate controls on vegetation phenological patterns in northern mid and high latitudes inferred from MODIS data. Global Change Biology, 10, 1133–1145. Zhang, X., Friedl, M.A., Schaaf, C.B., Strahler, A.H., & Schneider, A. (2004b). The footprint of urban climates on vegetation phenology. Geophysical Research Letters, 31(L12209). http://dx.doi.org/10.1029/2004GL020137. Zhang, P., Imhoff, M.L., Bounoua, L., & Wolfe, R.E. (2012). Exploring the influence of impervious surface density and shape on urban heat islands in the northeast United States using MODIS and Landsat. Canadian Journal of Remote Sensing, 38(4), 441–451. Zhang, P., Imhoff, M.L., Wolfe, R.E., & Bounoua, L. (2010). Characterizing urban heat islands of global settlements using MODIS and nighttime lights products. Canadian Journal of Remote Sensing, 36, 185–196. http://dx.doi.org/10.5589/m10-039. Ziska, L.H., Gebhard, D.E., Frenz, D.A., Faulkner, S., Singer, B.D., & Straka, J.G. (2003). Cities as harbingers of climate change: Common ragweed, urbanization, and public health. Journal of Allergy and Clinical Immunology, 111(2), 290–295. http://dx.doi.org/10. 1067/mai.2003.53.