Urban heat island impacts on plant phenology: intra-urban variability and response to land cover

Our study domain is Dane County (WI, USA), at the center of which is the city of Madison. Madison is a mid-sized urban area located in north-central United States (43°N, 89°W) with a population of 233 000 and an estimated urban agglomeration population of 402 000 (US Census Bureau 2010). The climate is humid continental, characterized by consistently below-freezing winter temperatures and warm summers with precipitation dominated by convection-based storms. Madison is surrounded by an agricultural landscape intermingled with deciduous forest and lakes (figure 1(a)).

Zoom In Zoom Out Reset image size

Beginning in March 2012, 151 HOBO U23 Pro v2 temperature/RH sensors (Onset Computer Corporation, Bourne MA) were installed on utility and streetlight poles in and around Madison. 135 sensors were installed in March 2012, with an additional six sensors installed in October 2012 and ten more in August 2013. All sensors were equipped with radiation shields, installed at a height of 3.5 m, and logged instantaneous temperature every 15 min. Sensors were classified as Urban, Rural, Park, Lake, or Wetland at the time of installation. Specifically, urban sensors were defined as within the municipal limits of one of the cities or towns within the Madison metropolitan area (including Madison, Monona, Middleton, Sun Prairie, Fitchburg, Verona, Waunakee, and Stoughton), while rural sensors were outside of these boundaries. Park sensors were within municipal boundaries of a city or town, and inside an officially designated park; parks range in size from 1.7 to 520 ha, with a mean (median) size of 101 ha (39.2 ha). A full description of the sensor network, which is one of the most spatially dense and extensive in existence, can be found in Schatz and Kucharik (2014).

A detailed description of our approach for estimating SOS/EOS based on temperature data is described in appendix S1. In short, we used a heating and cooling degree day (HDD and CDD) approach developed by Richardson et al (2006) for sugar maple, which is the Wisconsin state tree and an important component of Madison's urban canopy. The HDD/CDD approach is based on the concept of thermal time and estimates the start and end of growing seasons based on cumulative temperature above/below a physiologically meaningful threshold (Cannell and Smith 1983, Fu et al 2014a).

As this approach represents only a single species and was developed in New Hampshire, we also used a permutation-based approach to quantify uncertainty associated with our estimate, in which temperature thresholds and other necessary parameters are randomly selected from a distribution in order to provide a range of values, from which both a mean and uncertainty can be estimated (Serbin et al 2014, Zipper and Loheide 2014). We selected a uniform distribution from −0.6 °C to 4 °C as our temperature threshold range for SOS thresholds for HDD accumulation, as this interval encompasses a range of vegetation, from cold-tolerant to warm-season, and these temperatures have been demonstrated as effective phenological predictors previously (Schwartz and Marotz 1986, Richardson et al 2006). As less data exists for the autumn phenological triggers, we varied our EOS threshold for CDD accumulation over the same 4.4 °C range centered on the 20 °C value reported by Richardson et al (2006) as most effective for sugar maple.

For each sensor and each growing season, cumulative HDDs beginning January 1 and CDDs beginning August 15 were used to identify the temperature-based SOS/EOS (SOS_T and EOS_T, respectively), and GSL_T was calculated as EOS_T − SOS_T. The 2012 SOS_T is an exception to this, as the sensors were installed in late March and therefore temperature from a meteorological station at the University of Wisconsin-Madison Arboretum were adjusted for the UHI effect unique to each sensor and used prior to March 27. Differences between urban, park, and rural sensor classes were tested for statistical significance using a pairwise t-test with significance defined as p < 0.05.

A detailed description of the methodology for calculating phenology from remotely sensed data is presented in appendix S2. We collected all Landsat images of Madison WI (path 24, row 30) over the period 2003–2013 with < 40% cloud cover, for a total of 110 images (figure S1). Each image was atmospherically corrected using the LEDAPS toolchain (Masek et al 2006), cloud-masked using the FMASK utility (Zhu and Woodcock 2012), and clipped to the boundaries of Dane County. We used a double-logistic model (Zhang et al 2003, Fisher et al 2006) to calculate a composite remotely sensed SOS/EOS (SOS_R/EOS_R) based on five different vegetation indices and calculated GSL_R as EOS_R − SOS_R.

It is important to note that this timestacking method produces a representative or average SOS_R/EOS_R/GSL_R for the period 2003–2013 at each pixel, rather than SOS_R/EOS_R/GSL_R for each individual growing season (Fisher et al 2006, Elmore et al 2012, Melaas et al 2013). Due to infrequent Landsat overpasses (figure S1) and common summer cloud cover in the area, using Landsat imagery over a smaller time period such as 2012–2014 for direct comparison with temperature-based method is not possible. However, satellites with more frequent overpasses (e.g. MODIS, AVHRR) have spatial resolutions an order of magnitude coarser, which would make studying the fine-scale impacts of the UHI on phenology impossible. While data fusion from multiple sensors shows promise at obtaining phenological estimates at high spatial and temporal resolution (e.g. Liang et al 2014), fusion approaches currently struggle when pixels have a mix of land cover types, and therefore are not well-suited to urban analysis (Walker et al 2012, Zhang et al 2013, Klosterman et al 2014).

To address potential problems associated with comparing 2003–2013 remotely sensed imagery with 2012–2014 temperature data, all pixels containing land use change between 2001 and 2011 (Jin et al 2013) were masked from analysis. Furthermore, we separately fit and estimated phenology using data from the complete Landsat period of record (1982–2013) and carried out all subsequent analysis to determine whether choice of the Landsat period has an impact on results. Only results using the 2003–2013 Landsat imagery are presented here, as using 1982–2013 data did not substantially impact either results or interpretation (data not shown).

To compare remotely sensed and temperature-based phenology, we averaged SOS_R, EOS_R, and GSL_R within a 500 m buffer surrounding each sensor, as this distance has been recommended as a zone of influence for sensor-based studies (Oke 2006) and previous work found that 500 m best describes local climate variability in our study area (Schatz and Kucharik 2014). Due to the different time periods of remotely sensed (2003–2013) and temperature-based (2012–2014) methods, it is important to note that a 1:1 relationship of phenological indicators between the methods is not expected, as the temperature regime controlling the timing of SOS/EOS varies interannually. However, temperature-based results show a consistent year-to-year relationship between SOS_T/EOS_T and impervious cover (discussed in section 3.1), indicating that the spatial patterns of UHI impact on meteorological conditions conducive to plant growth are consistent from year-to-year. Therefore, we hypothesize that a positive correlation reflecting these patterns exists between SOS_T/EOS_T from individual growing seasons and long-term average SOS_R/EOS_R. Our analysis focuses on the correlation and slope of the relationship between remotely sensed and temperature-based methods rather than bias between the two methods, as bias likely reflects the unique characteristics of the 2012–2014 growing seasons compared to the 2003–2013 average (see Ault et al 2013).

As a first-order attempt to quantify the impact of land cover on differences between remotely sensed and temperature-based growing season, we define three difference metrics. For each difference metric, a positive value indicates that remotely sensed estimates correspond to a longer growing season than would be predicted from meteorological data (e.g. SOS_R < SOS_T, EOS_R > EOS_T, and GSL_R > GSL_T):

• (1)  
  SOS_Diff = SOS_{T_Norm −} SOS_{R_Norm}
• (2)  
  EOS_Diff = EOS_{R_Norm −} EOS_{T_Norm}
• (3)  
  GSL_Diff = GSL_{R_Norm −} GSL_{T_Norm}

To account for the different temperature-based growing seasons (2012–2014) within a single model, we calculated normalized values of SOS_T, EOS_T, and GSL_T (denoted SOS_{T_Norm}, EOS_{T_Norm}, and GSL_{T_Norm}) by subtracting the yearly mean value across all sensors from each individual sensor's value, thus centering our data at 0 for all years. Similarly, we subtracted the overall mean value from each sensor's value for the remotely sensed metrics SOS_R, EOS_R, and GSL_R to estimate SOS_{R_Norm}, EOS_{R_Norm}, and GSL_{R_Norm}. Using this technique, there are three difference metrics for each sensor (once for each growing season).

To account for potential impacts of land use change across the interval over which Landsat results were averaged, we eliminated from analysis all pixels classified as changing land cover between 2001 and 2011 in the National Land Cover Dataset (NLCD) (Jin et al 2013). Furthermore, we repeated our analysis using SOS_R and EOS_R derived from a different Landsat time window (1982–2013) to check for the impacts of time domain selection on mean phenology, which had no substantial difference in results or interpretation (only 2003–2013 results shown below).

We used two approaches to characterize the relative composition of impervious, tree, and grass cover at 30 m resolution (figures 1(b)–(d), both relying on spectral mixture analysis. Spectral mixture analysis is a widely used technique for estimating relative composition of different land cover classes by comparing the spectral signature of each pixel to spectra from user-defined reference endmembers (Small 2001, Wu and Murray 2003, Wu 2004, Buyantuyev et al 2007, Weng 2012). Inputs to both approaches were 30 cloud-free Landsat images collected between 2003 and 2013 with cropland and water masked (Gan et al 2014), and spectral mixture analysis was conducted in ENVI 5.0 (Exelis Visual Information Solutions, Boulder, Colorado) and constrained to a unit sum.

In the first approach (referred to as the SMA method), we performed a linear spectral mixture analysis, trained with areas of homogeneous cover manually selected in each of the three classes. In the second approach (referred to as the hybrid method), we use percent impervious cover from the NLCD (Xian et al 2011), a national-level impervious cover estimate widely used in the UHI literature (e.g. Georgescu et al 2012, 2013, Zhang et al 2014, Walker et al 2015, Schatz and Kucharik 2015, 2016). We then used spectral mixture analysis to differentiate between tree and grass cover for the remaining percent land cover in each pixel not accounted for by NLCD impervious estimates.

We assessed the accuracy of each method by comparing results to manually digitized 2013 Dane County imagery from the US Department of Agriculture's National Agricultural Imaging Program (NAIP). NAIP imagery was manually classified into 'birds-eye view' percent impervious, tree, and grass cover at 50 randomly selected 90 m × 90 m blocks, which were compared to 3 pixel-square blocks of Landsat-derived land cover (Wu and Murray 2003, Wu 2004, Gan et al 2014). Comparison between the two methods revealed only minor differences in results which did not alter overall interpretation, but a slightly better performance of the hybrid method; thus, results and analysis based on the hybrid method are reported in the main body text. Accuracy assessment results and a comparison of the two methods is presented in appendix S3.

The temperature-based growing season shows a strong influence of the UHI (figures 2 and 3). Median urban sensor SOS_T is advanced by 0.8, 0.9, and 2.0 d relative to the rural sensors in 2012, 2013, and 2014, respectively, and urban EOS_T delayed by 4.2, 3.2, and 3.3 d. These changes in SOS_T and EOS_T lead to a statistically significant median increase in GSL_T by 4.9–5.3 d in urban areas (figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The changes in SOS_T and EOS_T are related most strongly to the percent impervious land surface cover (%I), a proxy for the density of the built environment, in the area surrounding the sensor (figure 3). We averaged %I within 500 m of each sensor, as local climate variability in our study area is best described at this distance (Schatz and Kucharik 2014), but changes in buffer radius (100–1000 m) do not substantially alter results (data not shown). Across all years, R² for the %I–SOS_T/EOS_T relationships range from 0.51 to 0.68, indicating that physical urban density consistently explains over half of GSL_T variability. Uncertainty estimates generated using our permutation-based approach show similar trends as a function of impervious cover, indicating that variability in temperature thresholds or degree day accumulation intended to represent a range of potential urban species does not alter the relationship between the UHI and phenology.

Phenological shifts when considering rural and urban endmembers, which we define as the minimum (0.7%) and maximum (76.1%) observed impervious cover within a 500 m buffer of all sensors, were even more extreme than shifts in median SOS_T/EOS_T. The urban endmember SOS_T was advanced by 1.5–3.7 d and EOS_T delayed by 6.0–7.9 d relative to the rural endmember. This represents a potential UHI-driven extension of the growing season by 8.0–10.5 d, or 4.7%–6.2% across our three study years.

The magnitude of UHI-driven changes in phenology are strongly determined by the prevailing weather conditions during seasonal transitions. The 2012 and 2014 SOS_T provide a useful contrast to explore this dynamic. In 2012, much of the US including Madison WI experienced anomalously high late winter temperatures leading to a record-setting 'false spring' (Ault et al 2013). This was driven by a regional, synoptic-scale warm front which elevated temperatures above physiologically important thresholds more or less simultaneously across urban and rural parts of our study area, and therefore the small temperature differences caused by the UHI had relatively little effect. In contrast, 2014 experienced a relatively cool spring, which meant that small (1 °C to 2 °C) increases in temperature due to the UHI were able to contribute a larger proportion of the degree day requirements to trigger the onset of spring. Therefore, the magnitude of urban impacts on SOS_T were ∼2.4 times as strong in 2014 compared to 2012, as measured by the difference in slopes between the two growing seasons (− 0.020 d/% in 2012 and − 0.048 d/% in 2014). This indicates that studies reporting UHI-driven increases in GSL on the basis of a single growing season (e.g. Schmidlin 1989, Zhang et al 2004) may not be generalizable across time without considering prevailing weather conditions during spring and fall transition periods.

Park sensors tended to have intermediate values of SOS_T and EOS_T between the urban and rural sensors, indicating that they likely experience both a UHI effect elevating temperatures above the rural surroundings and a park cool island (PCI) effect reducing temperature relative to their urban surroundings (Taha et al 1991, Spronken-Smith and Oke 1999, Upmanis and Chen 1999, Spronken-Smith et al 2000, Yu and Hien 2006, Feyisa et al 2014). The PCI leads to a subdued expression of the UHI and median park GSL_T 1.7–3.2 d longer than rural sensors (p < 0.05 in 2013 and 2014) and 1.5–3.3 d shorter than urban GSL_T (p < 0.01 all years). This subdued phenological response to the UHI in parks may be particularly beneficial for native and migratory species which rely on vegetation phenological events for survival (e.g. budburst), as urban parks are often viewed as refugia for species threatened by the urbanization process (Dobrowski 2011, Stagoll et al 2012).

SOS_R, EOS_R, and GSL_R are shown in figures 4(a), (c), (e). Patterns of urban and peri-urban land surface phenology are consistent with those observed in other studies, with earlier SOS_R and later EOS_R observed close to the center of the city (Fisher et al 2006). Due to the elimination of agricultural lands surrounding Madison from analysis, phenological estimates are much less dense outside the urban margin. However, it is evident that GSL_R is shorter in rural landscape surrounding Madison, with several patches of longer GSL_R corresponding to small outlying cities and towns. Areas that are predominantly wetland (e.g. northeast and south of Madison) have the shortest observed GSL_R within our study domain, indicating that the growing season of naturally cool areas is well captured.

Zoom In Zoom Out Reset image size

At the county scale, GSL_R increases of ∼10–25 d are observed within ∼10 km from the city center (defined as the Wisconsin state capitol building; purple star in figure 4(e)) and maximum increases occur ∼2–7 km from the city center (figure 4(f); purple circles in figure 4(e)). These patterns exist due to the combined effects of the Madison lakes (suppressing GSL_R within 2 km of the city center) and the urban–rural transition (occurring ∼6.5–11 km from the city center).

We believe that GSL_R is suppressed within the 2 km closest to the city center due to a unique characteristic of Madison, which is that the densest parts of the city are on an isthmus between two lakes (figure 1). While these lakes visually dominate the landscape, previous work by Schatz and Kucharik (2014) has found that the lake's zone of influence on temperature is relatively small (100s of meters), and impacts of lakes on temperature decay exponentially with distance from the lakeshore. Therefore, the 2 km closest to the city center (equal to the half-length of the isthmus and shown as the inner purple circle in figure 4(e)) represents the area over which lakes contribute a larger proportional influence on temperature and phenology; however, previous work has shown the influence of the lakes on temperature is minor outside of this narrow zone on the isthmus (Schatz and Kucharik 2014).

Instead, the urban–rural transition represents the primary control over observed patterns in GSL_R, a pattern similar that shown by Fisher et al (2007) in Providence RI and supporting recent work indicating that urban form determines the spatial extext of urban impacts on the thermal environment (Yang et al 2016). The city of Madison extends ∼6.5 km south from the city center, and ∼8–11 km on the east and west sides of town; therefore, the radial interval from ∼6.5 to 11 km distance represents the urban–rural transition zone. This transition is evident in figure 4(f), where we can see a decrease in GSL_R over this approximate interval, and relatively static baseline rural GSL_R at distances greater than 11 km. The decrease in GSL_R over this urban–rural transition zone can be explained by two mechanisms. First, this transition is characterized by a decrease in impervious cover with increasing distance from the city center (figure 1). Based on the observed relationship between GSL_T and impervious cover (section 3.1), the meteorological period suitable for vegetation growth is longer in urban areas with denser impervious cover and a warmer temperature regime, leading to a lower GSL_R in the surrounding rural areas. Second, urban areas are often characterized by exotic, ornamental, or invasive vegetation (McKinney 2002, Niinemets and Peñuelas 2008), including non-native evergreen plants, which may be characterized by a longer growing season (Shustack et al 2009, Nemec et al 2011). Thus, we suggest that decreases in non-native vegetation moving outward from the city center contributes to the observed decrease in GSL_R over the urban–rural transition zone and reinforces previously described temperature effects.

While temperature-based and remotely sensed metrics show comparable patterns of longer growing seasons in urban Madison compared to the rural surroundings, when SOS/EOS are compared directly there is no statistically significant relationship between temperature-based and remotely sensed phenology at the individual sensor level (figures 4(b) and (d); R² < 0.06). As noted in section 2.4, we did not expect there to be a 1:1 relationship between SOS_T/EOS_T (which represents the individual 2012–2014 growing seasons) and SOS_R/EOS_R (which represents the mean of the 2003–2013 growing seasons); rather, we hypothesized that there would be a significant positive correlation between the two methods due to persistent interannual temperature effects on phenology described in section 3.1, with a bias between methods due to the different sampling intervals. To study the drivers of differences between temperature-based and remotely sensed metrics, we introduced three difference metrics (SOS_Diff, EOS_Diff, GSL_Diff) in section 2.4. For each difference metric, a positive value indicates that remotely sensed estimates correspond to a longer growing season than would be predicted from meteorological data (e.g. SOS_R < SOS_T, EOS_R > EOS_T, and GSL_R > GSL_T).

We find significant relationships between our different metrics and the land cover in the area surrounding each sensor (figure 5). Relationships between impervious cover (%I) and our difference metrics are best explained by a piecewise linear function, with a breakpoint determined as the %I which minimizes the sum of squared errors between model fit and difference metrics. Breakpoints are consistent across difference metrics at %I of 35% (SOS_Diff), 37% (EOS_Diff), and 36% (GSL_Diff) (figure 5). Beneath this breakpoint, increasing %I is associated with a longer realized growing season relative to meteorological potential at a rate of 6.5 d/10%, with approximately even contributions from changes to SOS_Diff and EOS_Diff. Above this breakpoint, increasing %I is associated with a shorter growing season at a rate of 4.1 d/10%, with ∼2 times stronger effects on EOS_Diff than SOS_Diff.

Zoom In Zoom Out Reset image size

We attribute this piecewise pattern to two drivers: vegetation composition and urban stress regimes. First, moving outward from the city center, land cover transitions between the following classes: (Zone 1) high-density urban areas dominated by impervious cover, to (Zone 2) old neighborhoods with a thick urban canopy dominated by tree cover, to (Zone 3) recently developed areas with sparse urban canopies and a stronger grass signal, to (Zone 4) rural areas predominantly characterized by agriculture and trees (figures 1(c)–(e)). Sensors with peak difference metrics (30% < %I < 40%) are primarily in low-density urban areas (Zone 3) surrounding Madison, as well as outlying towns and villages, where grass cover is highest, and the other sensors making up the increasing limb (%I < 36%) are in rural areas or the urban outskirts (Zone 4), where most remaining natural vegetation is trees. Therefore, the increasing limb of the %I relationships represents an increase in the relative proportion of cold-tolerant urban turfgrass within vegetated areas (%G), which is typically green from shortly after snowmelt in the spring until the first winter snowfall, moving from Zone 4 into Zone 3. This is supported by a significant positive correlation between %G and SOS_Diff (R² = 0.24, p < 0.001), as well as a weakly positive correlation with EOS_Diff (R² = 0.02, p < 0.05) (figure 5). Overall, GSL_Diff increases (i.e. the observable period of greenness grows longer relative to our temperature-based estimates) by 2 d for every 10% increase in %G, with ∼75% of that change occurring at the beginning of the growing season. The decreasing limb of the %I relationships, then, may also be partially explained by the decrease in grass cover associated with the transition from low-density urban areas surrounding Madison (Zone 3) to the higher-density areas closer to the city center (Zones 1 and 2) in which %G decreases.

Second, we suggest that the decreasing limb observed in %I relationships is may be due to increased water or pollutant stress in denser urban areas (Zones 1 and 2), which can lead to early onset of senescence and effectively decouples EOS_R from meteorological conditions (Gratani et al 2000, Honour et al 2009, Sjoman and Nielsen 2010). This is supported by the observation that the slope of the decreasing relationship between %I and EOS_Diff is approximately double the slope of the decreasing relationship between %I and SOS_Diff. While our analysis does not consider irrigation, variability in water available to plants as a result of urban irrigation (e.g. Pataki et al 2011, Bijoor et al 2012, Vico et al 2014) may further decouple EOS_R from EOS_T, e.g. by allowing irrigated vegetation to remain green during drought while non-irrigated vegetation may senesce.

While our analysis is conducted at a plant functional type level (e.g. grass versus trees), previous field- and plot-based studies have shown that phenological response to temperature varies at the species level (Bradley et al 1999, Chuine 2000, Primack et al 2004, Morin et al 2009, Vitasse et al 2009, Jochner et al 2013). A further complication is that non-native vegetation is common to urban areas and may be associated with earlier greening in the spring and a concomitantly longer GSL_R, as discussed in section 3.2 (McKinney 2002, Niinemets and Peñuelas 2008, Shustack et al 2009, Nemec et al 2011). If we assume that non-native vegetation common to urban areas has a longer growing season than native vegetation, this would cause an earlier SOS_T (downward shift of urban sensors in figure 4(b)) and a later EOS_T (upward shift of urban sensors in figure 4(d)), both of which would improve the correlation between temperature-based and remotely sensed phenological metrics.

While our temperature-based method accounts for potential variability in species-specific biological responses to thermal conditions by varying input parameters to our HDD/CDD equations (see appendix S1), this technique implicitly assumes a random and homogeneous distribution of all species by equally weighting all parameter combinations when calculating the mean SOS_T/EOS_T for each sensor. To account for spatial variability in species composition, the fitting parameters in table S1 could be 'tuned' to maximize agreement between SOS_T/SOS_R and EOS_T/EOS_R at each sensor (e.g. Fisher et al 2007); however, this approach would assume that temperature is the only factor contributing to spatial variability in phenology, whereas factors such as intra-species variability in response to photoperiod or environmental stressors may contribute to spatial variability in phenology (Gratani et al 2000, Saxe et al 2001, Schaber and Badeck 2003, Honour et al 2009, Caffarra et al 2011). Urban pixels are typically mixes of the built environment and vegetation (Gao et al 2006, Klosterman et al 2014, Lazzarini et al 2015, Walker et al 2015), leading to greater variability in both land cover and phenology than natural areas (Buyantuyev and Wu 2012).

These results provide the first synthesis of remotely sensed phenology with a dense ground-based urban temperature sensor network, and highlight a disagreement between the two methods which is associated with fine-scale variability in land cover. As previous work has reported substantial bias between remotely sensed and observed phenology, SOS_R/EOS_R are typically used to compare the relative timing of phenological events, rather than the exact dates (White et al 2009, Cong et al 2012, Xu et al 2014). However, here we show that remotely sensed observations of variability in land surface phenology cannot be used as a proxy for UHI intensity, as the vegetative response to meteorological conditions is dependent on the highly variable land cover at and surrounding a point. Even when focusing on the response of a single plant functional type to the UHI, such as forests in Fisher et al (2006) and Elmore et al (2012), we find that land cover (particularly %I) surrounding a pixel alters realized GSL and needs to be considered before conclusions can be drawn about the strength of the UHI. These results extend previous work done at larger spatial scales showing a disconnect between meteorological conditions and phenological response (Fisher et al 2007) and indicates that local processes, particularly land cover composition, must be considered as an important control over the vegetative response to changes in land cover in addition to changes in climate.