Interactions between urban vegetation and surface urban heat islands: a case study in the Boston metropolitan region

Our study region includes the Greater Boston metropolitan area (figure 1), encompassing ∼12 000 km² of residential and commercial land use in and around Boston, with surrounding rural areas composed of northern hardwood and oak/hickory forests inland and woody wetlands along the coast. Elevation in the study region ranges from 0 to 192 m above sea level and the climate is classified as humid continental. Precipitation is uniform throughout the year, but temperatures are highly seasonal with relatively warm summers and cold winters. In response, deciduous vegetation in the region has well-defined phenology that is primarily controlled by temperature (e.g., Richardson et al 2006, Friedl et al 2014).

Zoom In Zoom Out Reset image size

We used all available Landsat TM and ETM+ images from the United States Geological Survey with less than 90% cloud cover for two scenes that cover the Greater Boston metropolitan area (Path 12, Rows 30 and 31; see figure 1). Image digital numbers were converted to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Processing System algorithm (Masek et al 2006), and observations contaminated by clouds, cloud shadows, and snow and ice were identified and removed using the Fmask algorithm (Zhu and Woodcock 2012). The final data set included a total of 1175 TM and ETM+ images spanning the period from 1984 to 2013.

To map information related to green leaf phenology, we used the Landsat Phenology Algorithm (LPA; Melaas et al 2013) to estimate the long-term average and annual day of year associated with leaf emergence (start of spring; SOS) and autumn senescence (end of season; EOS) at 30 m spatial resolution. Because the LPA was originally designed for application over forested areas with relatively dense tree cover (i.e., with relatively high maximum summertime EVI), we incorporated two modifications to support detection of phenology in areas with lower tree cover (i.e., urban areas). First, instead of estimating separate spring and autumn logistic curves to model seasonal variation in Landsat EVI, we used smoothing splines to model the mean annual phenology at each location, and then normalized the resulting function at each pixel to vary between 0 and 1. Similar to Melaas et al (2013), the long-term average SOS and EOS at each pixel were estimated based on the date when the normalized spline reached 50% of the spring and autumn amplitude, respectively. Second, using a high-resolution canopy cover map of Boston (see below), we found that arbitrarily defined ranges for long-term mean winter background and peak summer EVI significantly underestimates deciduous tree coverage within the city. Alternatively, since trees are less moisture-limited than grasses and shrubs, they exhibit less interannual variability in summer EVI and, hence, greater agreement (i.e., higher correlation) with a long-term mean fitted cubic spline. Therefore, using the canopy cover map, we identified deciduous forest pixels in our study region where the Pearson correlation between observed EVI values and the fitted cubic spline was greater than 0.85 (figure 2). When the correlation was less than 0.85, the pixel did not consistently follow the relatively stable annual cycle of greenness of deciduous trees, and was assumed to contain mostly grasses and shrubs (e.g. lawns and golf courses) or ISA (see supplementary material).

Zoom In Zoom Out Reset image size

Remotely sensed SUHI intensity is widely mapped over various scales using land surface temperature (LST; Yuan and Bauer 2007, Peng et al 2011). For this work, we converted top-of-atmosphere brightness temperatures from the Landsat TM and ETM+ thermal band to LST for all images between May 1 and September 1 using a first-order thermal radiative transfer equation:

$$\begin{eqnarray*}&&{L}_{\lambda }=[\varepsilon {L}_{{\rm{T}}}+(1-\varepsilon ){L}_{{\rm{D}}}]\tau +{L}_{{\rm{U}}},\end{eqnarray*}$$

where $\tau$ is the atmospheric transmissivity, $\varepsilon \;$ is the land surface emissivity, ${L}_{\lambda }$ is the sensor-detected radiance, ${L}_{{\rm{T}}}$ is the target pixel surface radiance, and ${L}_{{\rm{D}}}$ and ${L}_{{\rm{U}}}$ are the downwelling and upwelling path radiances at the surface and top-of-atmosphere, respectively (Sobrino et al 2004). The radiance values (${L}_{\lambda }$ and ${L}_{{\rm{T}}})$ were converted to brightness temperatures (T_b) using sensor-specific Landsat calibration constants and the inverse Planck equation:

$$\begin{eqnarray*}&&{T}_{{\rm{b}}}=\displaystyle \frac{{k}_{2}}{\mathrm{ln}\left(\displaystyle \frac{{k}_{1}}{{L}_{{\rm{T}}}}+1\right)},\end{eqnarray*}$$

where k₁ = 666.09 W m⁻² sr⁻¹ μm⁻¹ and k₂ =  1282.71 K for Landsat 7 ETM+, and k₁ =  607.76 W m⁻² sr⁻¹ μm⁻¹ and k₂ = 1260.56 K for Landsat 5 TM (Barsi et al 2005). Atmospheric correction was performed using an algorithm based on the MODTRAN5 atmospheric radiative transfer model (Barsi et al 2005; http://atmcorr.gsfc.nasa.gov/), which uses atmospheric profiles of temperature, humidity, and pressure from the National Centers for Environmental Prediction (NCEP) to estimate values for $\tau ,$ ${L}_{{\rm{U}}},$ and ${L}_{{\rm{D}}}.$ We ran MODTRAN at 42° N, 71° W and at 15:16 GMT, which to corresponds to the approximate location of Boston and time of the Landsat overpass. Note that because the MODTRAN tool provided by Barsi et al (2005) includes NCEP data extending back to 2000, and the LST data used in this analysis includes imagery from 2000 to 2013.

To prescribe the emissivity at each pixel, we computed area-weighted emissivities based on three sub-pixel land cover constituents (vegetation, bare soil, and urban):

$$\begin{eqnarray*}&&\varepsilon ={f}_{{\rm{urb}}}{\varepsilon }_{{\rm{urb}}}+{f}_{{\rm{veg}}}{\varepsilon }_{{\rm{veg}}}+{f}_{{\rm{soil}}}{\varepsilon }_{{\rm{soil}}},\end{eqnarray*}$$

where the $f$ subscript denotes the sub-pixel fraction of urban, vegetation, and soil background within each pixel, $\varepsilon \;$ is the effective emissivity over the Landsat thermal bandpass, and where values for constituent emissivities (ε_urb = 0.9394, ε_veg = 0.9740, ε_soil =  0.9679) were based on published data (Baldridge et al 2009; see supplementary material). Values of fractional urban land cover, ${f}_{{\rm{urb}}},$ were estimated using the Massachusetts Office of Geographic Information's ISA layer (MassGIS 2007), which was generated using 50 cm near infrared orthoimagery and includes all constructed surfaces and man-made compacted soil (i.e., where no vegetation is present). To estimate the fraction of vegetative cover in each pixel, ${f}_{{\rm{veg}}},$ we used an approach similar to that used by Carlson and Ripley (1997) and Sobrino et al (2008) (who estimated ${f}_{{\rm{veg}}}$ as a function of NDVI), but used EVI because it provides greater dynamic range than the NDVI:

$$\begin{eqnarray*}&&{f}_{{\rm{veg}}}={\left[\displaystyle \frac{{\rm{EVI}}-{{\rm{EVI}}}_{{\rm{min}}}}{{{\rm{EVI}}}_{{\rm{max}}}-{{\rm{EVI}}}_{{\rm{min}}}}\right]}^{2},\end{eqnarray*}$$

where ${{\rm{EVI}}}_{{\rm{min}}}$ (=0.123) and ${{\rm{EVI}}}_{{\rm{max}}}$ (=0.660) were determined using a 1 m spatial resolution map of canopy cover map for Boston (Raciti et al 2014). The fraction of bare soil in each pixel, ${f}_{{\rm{soil}}},$ was then calculated as the residual of ${f}_{{\rm{urb}}}$ and ${f}_{{\rm{veg}}}$:

$$\begin{eqnarray*}&&{f}_{{\rm{soil}}}=1-{f}_{{\rm{veg}}}-{f}_{{\rm{urb}}}.\end{eqnarray*}$$

This formulation does not include water, which was excluded from the analysis. Further, in contrast to some studies that model sub-pixel cover as a combination of 'high albedo' and 'dark' urban materials (Small 2001, Small and Lu 2006), we assumed that areas located in the ${f}_{{\rm{urb}}}\;$ portion of each pixel were uniform in their emissivity.

Our approach for assessing SUHI impacts on green leaf phenology included two main elements. First, to provide a baseline characterization of SOS and EOS patterns across Boston's urban-to-rural gradient, we selected five 50 km transects originating at the Boston Common, a forested park in downtown Boston (figure 1). Two transects were located northeast and southeast of Boston (including some coastal areas), two were located inland to the northwest and southwest, and the fifth transect was located west-northwest of the city. Using these transects, we quantified how the timing of mean EOS and SOS vary as function of (1) distance from the city center and (2) average summertime LST, which we define here as the long-term mean LST based on all available clear-sky observations in May–August (2000–2013).

Second, we explored spatial patterns in SOS and EOS within the city of Boston, focusing on land use and vegetation controls. To characterize urban land use and vegetation cover, we estimated the sub-pixel composition of Landsat pixels within Boston's city limits for each of three surface types: impervious surfaces, tree canopy, and 'other' (grasses, wetlands, fields, bare soil etc). To estimate the fractions of ISA and tree canopy cover, we used the high spatial resolution maps of ISA and tree canopy described in section 2.3. The remaining surface area in each pixel was assigned to the 'other' class (figure 3). As part of this analysis, we also explored if and how two years with unusually cold and warm springs—2003 and 2010, respectively—affected patterns in SOS and EOS in Boston.

Zoom In Zoom Out Reset image size

Results from our analysis show that average growing season length (GSL; EOS–SOS) is approximately 18–22 days longer and mean growing season LST is 7 °C higher in urban areas relative to surrounding rural forests, and that spatial patterns in GSL and LST closely follow land use (figure 1). These relationships are consistent across Boston and in nearby cities (e.g., Lowell and Lawrence, located north of Boston in figure 1(a)). GSL is also significantly shorter for wetlands, most of which are located in coastal areas northeast and southeast of Boston.

LST and GSL patterns along the five transects in figure 1 show that relative to nearby rural locations, average SOS occurs about 10–12 days earlier (∼May 5) and EOS occurs about 8–10 days later (∼October 14) in locations that are within about 15 km of the city center (figure 4(a)). However, patterns in SOS and EOS also show substantial local variation (∼4–5 days) caused by heterogeneity in land cover and elevation. In contrast, the relationship between both EOS and SOS as a function of mean growing season LST is smoother (figure 4(b)), but has similar levels of local variation. On average, SOS occurs 1.8 days earlier and EOS occurs 2 days later for each 1 °C increase in mean growing season LST.

Zoom In Zoom Out Reset image size

Visual analysis reveals strong complementary gradients in ISA and tree canopy cover between northeastern and southwestern Boston, with pronounced local variation. Here we focus on fifteen neighborhoods and three green space areas (table 1, figure 3(d)). In the northeast, Charlestown (1), Downtown Boston (5), and East Boston (2) are heavily developed, but include areas with low stature vegetation (grass, shrubs) in parks along the coast and at Logan Airport in East Boston. Densely residential neighborhoods in Allston (6), Fenway (8), South Boston (3), and the South End (9) have moderate tree canopy cover and relatively high ISA. Lower density residential neighborhoods have much higher overall tree canopy cover, especially in Hyde Park (16), Jamaica Plain (12), and Roslindale (14).

Long-term mean SOS in Boston ranges from late April to late May (day of year 120–140), and long-term mean EOS ranges from early to late October (day of year 280–290) (figure 5). SOS occurs earlier in neighborhoods such as Jamaica Plain (11) and West Roxbury (15) that have moderate-to-high tree canopy cover mixed with ISA, and later in areas with very dense tree canopy cover such as Franklin Park (FP) and the Stony Brook Reservation (SB). Relative to SOS, spatial patterns in EOS are more variable and show less patch-structure. Spatial variability in both EOS and SOS is highest in neighborhoods such as Dorchester (11) that are the most heterogeneous with respect to tree canopy cover.

Zoom In Zoom Out Reset image size

SOS across Boston also shows distinct responses to interannual variation in temperatures. In 2003, for example, the average temperature for February, March, and April at Boston Logan International Airport was 2.9 °C below normal (4th coldest on record). In response, SOS anomalies in 2003 were generally positive (i.e., leaves emerged later than normal) by an average of 4.4 ± 3.4 days (±standard deviation; figure 5(c)). In contrast, temperatures during the same period of 2010 were 4.4 °C above normal (3rd warmest on record) and accumulated growing degree-days were nearly 150% above normal (results not shown). In response, SOS anomalies were uniformly negative (i.e., earlier leaf emergence) by an average of 10.7 ± 4.8 days (figure 5(d)). This effect was most pronounced in areas where average SOS tends to occur later such as the Stony Brook Reservation, where SOS was 15–20 days early in 2010 relative to long-term means.

Finally, to more fully explore the linkage between surface properties and spatial patterns in SOS and EOS, we tested the hypothesis that the amount of ISA in areas surrounding vegetation patches affects local SUHI intensity, and therefore influences the timing of SOS and EOS in urban vegetation patches. Results shown in figure 6(b) appear to confirm this hypothesis: SOS occurs 4.9 days later (t(d.f. = 25 580) =  75.4, p < 0.01) in areas surrounded by land cover with less than 20% ISA relative to locations in more developed areas. In autumn, EOS occurs 0.7 days later on average (t(d.f. = 42 761) = −12.1, p < 0.01) in areas where ISA is greater than 20%, but the relationship is relatively weak (see, figure 6(c)). In 2003 and 2010, pixels surrounded by areas with higher ISA leafed out about 1 day earlier and 2 days later (respectively) relative to normal (figures 6(e) and (f)), which suggests that the influence of ISA and SUHI on spatial patterns in SOS is weaker in years when spring temperatures are anomalous.

Zoom In Zoom Out Reset image size