Article Variations in Atmospheric CO2 Mixing Ratios across a Boston, MA Urban to Rural Gradient

Urban areas are directly or indirectly responsible for the majority of anthropogenic CO2 emissions. In this study, we characterize observed atmospheric CO2 mixing ratios and estimated CO2 fluxes at three sites across an urban-to-rural gradient in Boston, MA, USA. CO2 is a well-mixed greenhouse gas, but we found significant differences across this gradient in how, where, and when it was exchanged. Total anthropogenic emissions were estimated from an emissions inventory and ranged from 1.5 to 37.3 mg· C· ha −1 · yr −1 between rural Harvard Forest and urban Boston. Despite this large increase in anthropogenic emissions, the mean annual difference in atmospheric CO2 between sites was approximately 5% (20.6 ± 0.4 ppm). The influence of vegetation was also visible across the gradient. Green-up occurred near day of year 126, 136, and 141 in Boston, Worcester and Harvard Forest, respectively, highlighting differences in growing season length. In Boston, gross primary production—estimated by scaling productivity by canopy cover—was ~75% lower than at Harvard Forest, yet still constituted a significant local flux of 3.8 mg· C· ha −1 · yr −1 . In order to reduce greenhouse gas emissions, we must improve our understanding of the space-time variations and underlying drivers of urban carbon fluxes.


Introduction
The world's population has been rapidly shifting from rural and agrarian to urban areas, with the percent of global population living in cities increasing from 29.4% to 51.6% between 1950 and 2010 [1]. This urbanization trend is forecast to continue with models suggesting that by 2050 nearly 70% of the global population will live in urban areas [1]. While urban areas currently comprise less than 2% of global land area [2], their impact extends far beyond the city limits through environmental teleconnections [3] and demand for goods and services [4]. Urban areas are estimated to consume 67% of global energy and emit 71% of energy-related CO 2 emissions [5]. Despite the significant role these areas play in anthropogenic emissions, most research relating to atmospheric CO 2 dynamics has avoided areas close to or heavily influenced by cities [6,7]. Efforts to quantify terrestrial carbon exchange have instead focused on areas dominated by biogenic fluxes and homogenous land use patterns such as forest and agriculture [8]. By contrast, urban areas are often comprised of heterogeneous land cover and complex topography, which complicate measurements and source attribution of both CO 2 fluxes and mixing ratios.
A range of environmental gradients has been observed between urban and adjacent rural locations. For example, urban heat islands, where temperatures can be several degrees higher than adjacent rural areas, develop due to reductions in latent heat fluxes and surface albedo changes associated with paving, among other reasons [9,10]. Urban canyons created by tightly spaced buildings and roadways can change airflow patterns and increase downwelling longwave radiation by reducing the sky view factor. This in turn raises temperatures. Importantly, these increases in temperature have also been shown to extend the growing season [11,12] and likely also affect biogenic carbon exchange in cities. Differences in hydrology, floral and faunal species diversity, soil nitrogen and carbon stocks, and concentrations of atmospheric pollutants have also been observed along urbanization gradients [13], although not always according to expectations.
While some of the environmental gradients associated with urbanization have been better defined [14], the influence of these variables on atmospheric CO 2 mixing ratios has just recently begun to be assessed. For example, CO 2 mixing ratios have been found to be higher in urban centers compared to adjacent rural locations in Phoenix [15], Salt Lake City [16], and Baltimore [17]-a phenomenon known as an "urban CO 2 dome". These higher mixing ratios are due in part to local traffic emissions, as seen in Helsinki [18], Mexico City [19], and Basel [20], but may also be effected by residential, commercial, and industrial emissions. Unique patterns of CO 2 across urbanization gradients have also been demonstrated in Melbourne [21], Phoenix [22], and Rome [23], suggesting an association between urban land uses, urban density, and observed CO 2 [24]. There have also been attempts to map emissions at finer spatial scales in urban areas using mass flux measurements of carbon dioxide, among other data sources [25][26][27], but these results can be very difficult to interpret due to complex urban micrometeorology [28].
Urban areas' influence on atmospheric CO 2 is often framed in terms of anthropogenic emissions; however, biomass and biogenic CO 2 flux in urban areas can approach that of nearby forest-dominated areas [29][30][31]. In remote sensing products such as MODIS NPP, urban areas are masked out and assumed to have little productivity, but the biomass present in these areas suggests that biogenic fluxes are also important. Ecological processes in human-dominated ecosystems such as urban areas are expected to differ from adjacent, predominantly rural locations [32][33][34], but these differences in ecosystem function are poorly understood.
While our knowledge of carbon emissions and biogenic carbon exchange in urban areas is limited [31,35], local policies for climate action plans, emissions reductions, and urban greening are continually being developed (e.g., US Mayors Climate Protection Agreement and the California Global Warming Solutions Act). Biogenic carbon exchange estimates in urban areas are poorly constrained and represent a serious impediment to sustainable urban planning [36]. It is difficult to quantify the carbon exchange impacts of local greening initiatives such as Million Trees NYC and Grow Boston Greener, which have significant financial costs associated with them. Policymakers require better spatially and temporally resolved estimates of both anthropogenic emissions and biogenic exchange to assure that local climate mitigation actions are cost effective and CO 2 reductions are being actualized in the atmosphere.
In this study, we report atmospheric results from a new interdisciplinary research effort focused on (1) better characterizing atmospheric CO 2 mixing ratios across the urban-to-rural gradient near Boston, MA and (2) associating atmospheric CO 2 with changing CO 2 fluxes and land cover. We measured CO 2 mixing ratios at Harvard Forest in Petersham, MA (a forested area), Worcester, MA (urbanized to the east and forested to the west), and Boston, MA (urbanized) during 2011 in order to capture the heterogeneity of the urban gradient. From these observations, diurnal and seasonal patterns are examined. These temporal patterns are then compared to estimates of biogenic and anthropogenic CO 2 fluxes and land cover at each study site. We use remote sensing to investigate the potential implications of the urban heat island effect on vegetation phenology and atmospheric CO 2 exchange across the urbanization gradient. Finally, we explore the relationship of land cover to atmospheric CO 2 concentrations and determine how variables such as impervious surface area influence patterns of observed atmospheric CO 2 .

Results and Discussion
CO 2 is a well-mixed gas in the atmosphere; its spatial and temporal variations reflect a combination of anthropogenic emissions, exchange with the biosphere, atmospheric transport, and boundary layer dynamics. On the basis of atmospheric mixing patterns alone, CO 2 is expected to build up during the nighttime hours due to atmospheric stratification and decrease in the morning with the break-up of the nocturnal boundary layer [7,37]. However, the mixing ratio also reflects biogenic uptake, ecosystem respiration, and anthropogenic emissions including human respiration, each of which has different diurnal, seasonal, and spatial patterns. The biogenic signal tends to draw down daytime CO 2 during the summer growing season when photosynthesis is active, resulting in lower overall CO 2 mixing ratios. CO 2 mixing ratios are higher during the winter months when respiration dominates in the biosphere and heating-related emissions are highest.
These diurnal and seasonal trends in CO 2 mixing ratios were evident at each of our measurement sites, which spanned a gradient of urbanization intensities. All three sites had a larger seasonal amplitude in CO 2 than the global background measurements from Mauna Loa, HI and the measurements from Niwot Ridge [38], Colorado, which is a site at a similar latitude to our study area and within the free troposphere. These differences highlight broad scale patterns such as the increasing strength of seasonality with distance from the equator and the influence of local to regional uptake and release processes. In all cases, seasonal maxima and minima occurred during winter and summer months, respectively (Figure 1(a)). Total anthropogenic emissions estimates [39] also have a strong seasonal signal due to residential heating demand, and were roughly 4 and 24 times higher in Boston than similar estimates in Worcester and Harvard Forest, respectively (Figure 1(b)).

Trends in Observed CO 2
The spatial and temporal variability in carbon fluxes hinders simple characterization of the primary determinants of local CO 2 observations, especially in urban areas where the land cover is heterogeneous and topography is complex. Despite these challenges, we observed variations in CO 2 mixing ratios across Boston's urbanization gradient that were consistent with vegetation and urbanization patterns at each site ( Figure 2). For the Worcester site, which is the midpoint in our urbanization gradient, we separately characterized the results for air originating from the urbanized area to the east (0 to 180 degrees) and the rural area to the west (180 to 360 degrees) of the site. Henceforth, we will refer to these urban and rural sectors as East Worcester and West Worcester, respectively. Canopy percentage and biogenic fluxes (both  and ) were estimated within the 1 km radius around each tower (red circles). In Worcester, statistics were split into easterly and westerly sectors that represent half the areal coverage and are delineated by the dashed line. GPP = Gross primary production, E.R. = Ecosystem respiration, Human = Human CO 2 respiration, Mob. = Mobile source emissions, Res. = Residential emissions, and Other = All other fossil fuel emissions. All fluxes are in mg· C· ha −1 · yr −1 . All atmospheric CO 2 measurements are time weighted annual means with bootstrapped 95% confidence intervals.
In 2011, mean CO 2 mixing ratios in Boston were 8.8 ppm higher than air originating from East Worcester, 15.5 ppm higher than air originating from West Worcester, and 20.6 ppm higher than observations at Harvard Forest. These observations were consistent with the patterns in local anthropogenic and biogenic fluxes. Across all sites these differences amounted to a roughly 5% difference in mean annual CO 2 mixing ratios, despite the combination of a large biotic imprint on atmospheric CO 2 in the rural areas and large anthropogenic emissions in the urban areas. The 2011 annual mean observed CO 2 mixing ratios in Boston, East Worcester, West Worcester, and Harvard Forest were 393.4 ± 0.15, 398.5 ± 0.23, 405.2 ± 0.45, and 414.0 ± 0.21 ppm, respectively ( Figure 2). The trends in CO 2 mixing ratios at all sites showed seasonal shifts with winter enhancement-associated with heating related emissions and ecosystem respiration-and summer draw-down, coinciding with enhanced ecosystem productivity and reduced anthropogenic emissions.
Anthropogenic emissions for all sectors decreased significantly from Boston to Harvard Forest ( Figure 2). Total annual estimated fossil fuel emissions for Boston (excluding the area covered by water), East and West Worcester, and Harvard Forest were 34.7, 5.9, 1.97, and 1.53 Mg C· ha −1 · yr −1 , respectively. The composition of anthropogenic emissions also changed across Boston's urbanization gradient: emissions from other sources (such as industrial and commercial) decreased as a percentage of total emissions as urbanization decreased. There were also large seasonal differences in anthropogenic emissions: the ~47% increase in total emissions between summer and winter at all sites ( Figure 1) was driven by the ~480% increase in residential emissions between these seasons.
Patterns in the estimated biogenic fluxes showed the opposite trend; gross primary productivity (photosynthesis) and ecosystem respiration increased from urban Boston to rural Harvard Forest ( Figure 2). This increase in biogenic fluxes was associated with the increase in forest canopy from east to west across the region. Biogenic fluxes dominated carbon exchange processes at Harvard Forest and West Worcester and reflect the largely undeveloped, forested character of these areas. Tree canopy cover was 27% within a 1 km radius of the Boston tower: this is consistent with the 29% average overall canopy for the City of Boston [40]. Gross primary production and ecosystem respiration each constitute a substantial portion of total CO 2 fluxes in Boston, suggesting considerable biotic influence even within dense urban areas.
Differences between the human and vegetation-dominated environments across Boston's urbanization gradient were reflected in the annual standard deviations of CO 2 mixing ratios in Boston (17.8 ppm), East Worcester (21.5 ppm), West Worcester (15.9 ppm) and Harvard Forest (14.0 ppm). Higher overall and diurnal variability in Boston and East Worcester was likely due to proximate anthropogenic emissions, such as local traffic, combined with higher air entrainment from surrounding buildings. This variability was also exhibited in the hourly average CO 2 mixing ratios and the corresponding seasonal trends ( Figure 3). In Harvard Forest, total carbon emissions were relatively low in winter due to low biogenic and anthropogenic fluxes, resulting in CO 2 mixing ratios that remained relatively constant over time. Moving towards Boston, CO 2 mixing ratios quickly became more sinusoidal and reflected greater levels of anthropogenic emissions, which were quite high during winter months. The variability in the trend of seasonal CO 2 is supported by the changes in the heteroscedasticity and skewness of the frequency distributions of CO 2 observations at these three sites ( Figure 3). For example, the data distributions from Boston and East Worcester exhibited a strong positive skewness of +2.8 and +1.1, respectively. On the other hand, hourly CO 2 mixing ratios in West Worcester and Harvard Forest had a slight positive (+0.6) and negative skewness (−0.8), respectively, and exhibited a much lower variance. Without large proximate anthropogenic emissions at these two sites, mixing ratios rarely exceeded 450 ppm. The negative skew at Harvard Forest likely resulted from strong photosynthetic activity.
While other studies have observed urban CO 2 mixing ratios well above background levels, the magnitudes of these CO 2 domes varied greatly by location [15,21,41]. For example, mean peak city-center mixing ratios in Phoenix, AZ were 28%-76% higher than local background values, although this finding was likely influenced by highly stable atmospheric conditions resulting from local wintertime atmospheric inversion. In Portland, OR and Melbourne, Australia mean CO 2 mixing ratios at more developed sites were as much as 6 and 12 ppm greater, respectively, than those in corresponding lesser-developed locations. The strength of urban CO 2 domes, including Boston's, is sensitive to local meteorological conditions, emissions, biogenic processes, and the height of the gas analyzer above the surface. These local influences complicate simple generalization or extrapolation of urban carbon domes.
Observed diurnal patterns in CO 2 mixing ratios across Boston's urbanization gradient exhibited predictable behavior associated with stratification of the atmosphere, but also showed marked differences as urbanization increased ( Figure 4). The diurnal patterns at Boston and Worcester broadly showed a daily maximum occurring between 4:00 am and 7:00 am, followed by a rapid decrease occurring with sunrise and the associated break-up of the nocturnal boundary layer [20,37]. The daily minimum in CO 2 occurred in the early afternoon hours as atmospheric mixing and photosynthesis were maximized. At Harvard Forest, this same diurnal pattern occurred during the summer when photosynthesis and respiration were both large, but was absent during the winter months when CO 2 hovered around 400 ppm, reflecting the low local anthropogenic emissions (Figure 2), minimal photosynthesis, and reduced ecosystem respiration due to low temperatures [42]. We observed the largest diurnal variability in CO 2 during the summer months with maximum diurnal amplitudes of 29.2, 31.6, 31.1, and 29.0 ppm at Boston, East Worcester, West Worcester, and Harvard Forest, respectively ( Figure 4). These results vary slightly from observations across Portland, OR's urbanization gradient during summer and fall where amplitudes were higher in rural (33 ppm) and suburban (29.5 ppm) areas compared to the downtown core (25 ppm) [41]. Differences in both the absolute magnitude in CO 2 mixing ratios and their relationship with urban development between Boston and Portland's urbanization gradients reflect local meteorology, emissions, and the influence of deciduous versus evergreen vegetation exchange dynamics. For example, the Portland area has a greater number of conifers and a more temperate climate than Boston, which could result in biogenic fluxes that are greater in magnitude and driven more by moisture availability.
Diurnal patterns in mobile and total emissions in Boston and Worcester (east and west sectors combined) reflected human activity with overall emissions increasing around 7 am and remaining high through 8 pm ( Figure 5). Mid-day weekday CO 2 mixing ratios in Boston and Worcester were 5.1 and 2.3 ppm greater than on weekends, respectively. There was no statistically significant weekend effect observed at Harvard Forest. When integrated across the day, Vulcan mobile source emission estimates were 42.7% and 58.7% higher during the weekday compared to weekends in Boston and Worcester, respectively, which is consistent with elevated CO 2 mixing ratios observed during weekdays at each site. Observational studies in Portland and Phoenix showed weekday/weekend differences as high as 4.0 and 14.4 ppm, respectively [15,41], while a study from suburban Baltimore showed no significant weekend difference [17]. Weekend effects reflect the importance of local commuting patterns on observations of CO 2 mixing ratios. Despite being a relatively well-mixed gas, the imprint of human and biogenic activity can be seen in both the short-term and long-term signals of CO 2 across Boston's urbanization gradient. Moreover, many of the changes in CO 2 mixing ratios across the gradient were caused by alteration of land cover and the concomitant changes in vegetated fraction and anthropogenic emissions, as seen in Figure 2. These data suggest significant, direct alteration of CO 2 mixing ratios due to urban land cover change and associated anthropogenic activities.

Enhanced Vegetation Index (EVI) Time Series and Phenology Timing
While emissions variations are clearly associated with urban areas, less direct effects of urbanization can also influence CO 2 fluxes and observed CO 2 mixing ratios. The urban heat island effect, in particular, may alter the balance of respiration and photosynthesis in urban areas relative to nearby rural counterparts [11]. Urban heat islands may also alter seasonal anthropogenic emissions due to changes in heating and cooling degree days. Temperature gradients are frequently observed between urban and rural areas [9,43]. We observed mean summer temperatures of 21.7, 21.1, and 20.6 °C in Boston, Worcester, and Harvard Forest, respectively. Increases in temperature across the gradient were likely related to elevation differences, increased incoming longwave radiation (due to a reduced sky view factor and by the presence of atmospheric pollution), decreased latent heat fluxes, increased building and road storage heat flux, and anthropogenic heat emissions. [9,10,44].
Higher temperatures associated with urbanization can result in altered vegetation phenology [11,45]. For example, urban green-up (defined here as 25% leaf emergence) and brown-down (defined here as 90% leaf senescence) tend to occur earlier and later, respectively, than in nearby rural and suburban areas [12]. To examine trends in phenology across the gradient, we used the Enhanced Vegetation Index (EVI), which measures surface greenness. Based on the absolute EVI time series, green-up in 2011 occurred on approximate day of year (DOY) 126, 136, and 141 in Boston, Worcester and Harvard Forest, respectively ( Figure 6). Brown-down in 2011 occurred on approximate DOY 304, 284, and 288, respectively. The total growing season length difference between Boston and Harvard Forest in 2011 was 31 days, a potential 20% lengthening in the period for biogenic carbon uptake.
While the differential impacts of earlier green-up and later brown-down are still being quantified, for each one-day increase in growing season length, net ecosystem carbon uptake has been found to increase by 4.3 g· C· m −2 · day −1 across a range of temperate deciduous forests [46]. Assuming similar productivity per unit canopy cover and a 31 day phenologic change, the extended urban growing season could potentially increase net biogenic carbon sequestration in Boston (27% canopy cover) by as much as 0.36 mg· C· ha −1 · yr −1 , a 50% increase in net biogenic exchange. However, abiotic growing conditions that affect ecosystem productivity, such as soil moisture, nitrogen availability, and atmospheric ozone, differ significantly between urban and rural areas [13]. As a result, scaling ecosystem productivity by canopy cover should only be considered a first order estimate of the effect of longer growing seasons on carbon uptake. Further, while it is clear that the urban heat island effect significantly alters air and soil temperatures and growing season length in the region, it is difficult to determine what fraction of the lengthened growing season in Boston, relative to Harvard Forest, is due to heat island effects versus local climate and topographic differences between the sites.
The timing of green-up and brown-down occurred at different points along the seasonal CO 2 trend (Figure 3). In Boston, CO 2 mixing ratios began to decline in early February, well before the onset of photosynthesis in early May. This is likely due to a combination of increased vertical mixing, changes in background CO 2 mixing ratios, and the 50% reduction in residential emissions from January through March [39] due to the typical early year decrease in heating degree days (HDD): during 2011 in Boston, there were 1156 HDD in January, 961 in February, and 804 in March [47]. As a result, green-up occurred well after the timing of peak CO 2 mixing ratios in Boston. Conversely, CO 2 mixing ratios at all sites began to rise between late June and mid-August, proceeding brown-down by as many as 104 days in Boston's case. The late summer increase in CO 2 was likely not due to changes in residential emissions, but rather decreases in canopy photosynthetic efficiency (associated with foliar aging and decreasing insolation [42]), increasing ecosystem respiration, and changes to background CO 2 mixing ratios. Others have found that rates of gross primary production typically begin to decline in early July [48], consistent with the patterns observed at all sites. In contrast to the patterns in Boston, green-up and brown-down occurred closer to maxima and minima mixing ratios in Harvard Forest. Mixing ratios at Harvard Forest began to decline in mid-May, 25 days before green-up. Biogenic fluxes dominate this rural site, which likely reflects the influence of the coniferous trees in the canopy: conifers begin photosynthesizing as soon as daily mean temperatures are consistently above freezing [42]. The high values of EVI observed in Boston were influenced by urban and suburban lawns, which typically sequester less carbon than forests [49].

Land Cover and Relationship to CO 2
The differences we observed between our three study sites suggest that the land cover around a measurement location can influence observed mixing ratios of CO 2 . For example, we found that observations near a densely populated area with high traffic emissions (Boston) exhibited higher atmospheric mixing ratios than a site surrounded by forest (Harvard Forest). However, measurements in heterogeneous urban or suburban areas reflect a mosaic of land covers and depend on local meteorological conditions.
To better quantify the influence of land cover on atmospheric CO 2 , we conducted a more detailed analysis of the Worcester site, which has both large tracts of forests and urban development nearby (Figure 7). To the west of the Worcester site, forests dominate the land cover. To the east, residential, commercial/industrial, and other developed land uses dominate. Downtown Worcester-including Interstates 90, 190, and 290-is located between 40 and 170 degrees relative to the Worcester tower.
Total impervious surface area (ISA), which includes buildings, roads, and compacted man made soils, and mean EVI reflect these land covers (Figure 7).

Figure 7.
A basemap with EVI of the two sectors at Worcester tower site, derived from Zhu and Woodcock, 2012 [50]. Worcester was chosen for this analysis due to its proximity to large tracts of forest to the west and to Worcester's urban core, visible to the southeast. The 1 km and 5 km radii test areas were used to estimate biogenic emissions and correlate surrounding land cover to CO 2 observations, respectively. CO 2 mixing ratios measured from the primarily forested sector to the west (between 180 and 360 degrees) were on average 6.7 ± 1.8 ppm lower than CO 2 mixing ratios observed from the urban sector to the east (between 0 and 180 degrees) ( Figure 2). This provides further evidence that biogenic and anthropogenic fluxes associated with different land cover types likely influence observed atmospheric CO 2 mixing ratios. For example, West Worcester exhibited much higher tree canopy cover (65%) and corresponding estimated biogenic fluxes (−9.1 and +7.48 mg· C· ha −1 · yr −1 for GPP and E.R., respectively) than East Worcester (46% canopy, −6.44 and +5.29 mg· C· ha −1 · yr −1 for GPP and E.R., respectively). Anthropogenic emissions for East Worcester (11.16 mg· C· ha −1 · yr −1 ) were relatively high compared to West Worcester and Harvard Forest, suggesting that lower levels of canopy cover can serve as a proxy for human activity and associated anthropogenic emissions [51].
There were also seasonal differences exhibited between the sectors. The mean summer CO 2 mixing ratio for West Worcester was 389.4 ppm, compared to 395.2 ppm for East Worcester. These lower values in the west sector suggest increased photosynthesis, which is reflected in the higher EVI (0.56) and lower ISA (13.6%) exhibited by the west sector. Despite higher ISA (44.0%) and lower EVI (0. 41) in East Worcester, the estimated ecosystem respiration comprised the largest single source of CO 2 to the atmosphere (3.1 mg· C· ha −1 · yr −1 ), underscoring the importance of vegetation cover on local carbon fluxes, even within developed areas. In fall, the mean CO 2 mixing ratios for West and East Worcester were 404.0 and 416.1 ppm, and in winter, they were 407.0 and 419.6 ppm, respectively. The larger difference in mixing ratios between the two sectors in the fall and winter could be a result of proportionally greater increases in residential anthropogenic emissions in East Worcester relative to less-developed West Worcester during the cooler months.
While the land cover in Boston and Harvard forest was more uniformly urban and rural, respectively, we also parsed the CO 2 data by easterly and westerly wind sectors (Supplementary Material). In contrast to the Worcester results, the Boston and Harvard sector analyses in did not yield statistically significant differences, suggesting that the Worcester results were not driven by synoptic scale pollution patterns.
Despite these strong associations between land cover, CO 2 flux, and observed mixing ratios at the sector scale (east versus west), a more detailed spatial analysis, where we divided the source area into 45 degree wind sectors, yielded inconclusive results. There were no statistically significant correlations between ISA, EVI, or any of 8 land use classes and CO 2 mixing ratios for these more spatially resolved sectors. These results highlight several of the challenges associated with assessing the influence of heterogeneous land cover on CO 2 mixing ratios. Mean EVI was not necessarily lower in areas with higher human population and emissions: for example, the residential areas found in East Worcester (Figure 7) had both a high mean EVI and population density. We attribute these high EVI values in residential areas in part to grass, which increases mean EVI, but does not sequester as much carbon in biomass as forests [49]. This likely contributed to the inconclusive results in our analysis, further highlighting the challenges in using remotely-sensed measures of greenness as a proxy for biogenic fluxes in areas dominated by a mix of trees and grasses. Impervious surface area percentage is also a problematic proxy for CO 2 emissions. ISA can underestimate carbon emissions from point sources, which are small in area, and from major transportation arterials, which are relatively narrow and are often surrounded by vegetation to provide a sound buffer in urban areas. Furthermore, attribution of emissions to local energy usage is very challenging since energy demand is often spatially separated from power generation.
While we generated mixed results using this directional analysis around the Worcester site, there appeared to be an association between land cover and CO 2 mixing ratios across Boston's urbanization gradient, as demonstrated by Figure 2. As ISA and anthropogenic emissions increased from Harvard Forest to Boston, so did CO 2 mixing ratios.  .s.l.). The tower is located within a heavily developed urban area ~0.16 km from Interstate 90 and ~4 km from the downtown high-rise buildings. The BU Law building (~60 m a.s.l.) and the Warren Towers (~45 m a.s.l.) are approximately 70 and 180 m from the test site, respectively. Boston, with a population of approximately 600,000 people, is characterized by high-density urban development with parks interspersed. There are three power plants (all utilizing natural gas or oil) and one large regional airport within 16 km of the tower site [52].

Figure 8.
Land cover across our eastern Massachusetts study area, derived from Mass GIS data layers [53]. Impervious surface fraction is between 0 (no constructed impervious surfaces) and 1 (completely covered by constructed impervious surface). The base of the tower is approximately 349.4 masl. This area is characterized by low human population and by a mixed broadleaf forest dominated by red oak (Quercus rubra) and red maple (Acer rubrum) with moderately high biomass (115 mg· C· ha −1 ) [56]. The most proximate source of anthropogenic CO 2 emissions is Route 2, located ~5 km north of the tower. There are no large point sources of CO 2 within the Harvard Forest study area [52].
For this study, we examined meteorological observations and CO 2 mixing ratios at all three tower locations for 2011. We considered the biogenically dominated CO 2 mixing ratios from Harvard Forest to be the background signature for Eastern Massachusetts. Worcester and Boston represented points of increasing anthropogenic influence along Boston's urbanization gradient. The Boston, Worcester, and Harvard Forest areas had population densities of approximately 4,900, 2,800, and 9 people/km 2 , respectively [57]. Viewed together, these three sites constituted both an urban to rural gradient and an anthropogenic to biogenic carbon flux gradient.

Instrumentation and Data Analysis
In Boston, CO 2 and H 2 O mixing ratios were measured at 1 Hz using a Picarro 2301 cavity ring down spectrometer (Picarro, Sunnyvale, CA, USA). CO 2 measurements were corrected for water vapor content [58]. To ensure data quality and correct any possible sensor drift, periodic calibrations were performed every four months using three known reference CO 2 standards between 350 and 460 ppm (traceable to NOAA/CMDL). Instrument coefficients were adjusted to correct for minor drift over time (less than 0.2 ppm between calibrations). A Campbell Scientific (Logan, UT, USA) CSAT3 sonic anemometer was used to measure wind velocities and wind direction. Temperature and relative humidity were measured with a Vaisala HMP45C probe (Helsinki, Finland). The system operated nearly continuously from January through December: 2011 data completeness was 93%.
In Worcester, CO 2 mixing ratios were measured at 1 Hz with a closed path LI-6262 (LICOR, Lincoln, NE, USA) infrared gas analyzer. The system was automatically calibrated every six hours using three reference gases between 340 ppm and 460 ppm (traceable to NOAA/CMDL standards). Wind speed and direction were measured with a Met One 034B Windset 2 wind instrument. All data were recorded using a Campbell Scientific (Logan, UT, USA) CR1000 datalogger. The system operated nearly continuously from installation in late March through December: overall 2011 completeness was 70%. CO 2 at Harvard Forest was measured by a LI-6262 (LICOR, Lincoln, NE, USA) infrared gas analyzer. Automated calibrations of the sensor were run at least twice daily to account for instrumentation drift; at least two reference gases between 340 and 460 ppm (traceable to NOAA/CMDL) were used during these calibrations. Wind speed and direction above the canopy were measured using an Applied Technologies (Longmont, CO, USA) SATI/3K 3-D sonic anemometer. Air temperature was derived from sonic anemometer's speed of sound measurement by accounting for influence of water vapor. Relative humidity was measured by Vaisala (Vantaa, Finland) HMP45 probe in an aspirated radiation shield. Ambient atmospheric pressure was measured by a MKS Instrument (Andover, MA, USA) absolute manometer. Data at the EMS tower were digitized and recorded using a custom data acquisition and control system. The Harvard Forest EMS tower operated from January through December in 2011, but experienced several interruptions in the CO 2 data: 2011 data completeness was 65%.
Since instrument height can significantly alter observed CO 2 concentrations due to plume dispersion [28,59], gas analyzers were placed at approximately the same height above ground level at each location. While this ensured that all sites observed emissions from similarly sized source areas, each instrument was likely located at a different vertical position within its respective boundary layer. For example, observations from urban Boston and Worcester were likely affected by building topography and the corresponding changes in micrometeorology and atmospheric mixing [60,61]. Methodological challenges (e.g., urban tall tower construction) in elevating instruments above the roughness layer combined with differences in boundary layer height at each location prevented similar instrument positioning within each plume/boundary layer. However, given their locations on roofs of similar heights, the instruments in urban Boston and Worcester were assumed to be within a similar urban canopy layer, which is supported by the similar mean wind speed observed at these two sites.
On the same roof as the Worcester tower, there were also a series of large mechanical units. While we confirmed that these units were not venting fossil fuel exhaust, we also conducted wind rose analyses to verify that CO 2 -rich air from the interior of the building was not being vented (Supplementary Material).
The R software package, version 2.15.2, was used for all statistical analyses and for data pre-and post-processing [62]. Half hourly block averages were calculated from the quality-controlled data from each site. Given the non-random distribution of data gaps and strong seasonality in the CO 2 signal at all of the sites, annual means were calculated using time weighting such that the monthly means were first calculated then averaged to annual scales for 2011. Due to the non-normal data distributions, a bootstrapping method was used to determine the 95% confidence intervals (CI) of different ecosystem variables for the nine sample classes [63]. Unless noted otherwise, all parenthetically reported values are 95% CI. The seasonal trends in CO 2 were calculated using locally weighted regression (LOESS) with a 0.5 span [64].

Anthropogenic and Biogenic Carbon Fluxes
With the exception of human respiration, anthropogenic emission estimates for this analysis were based on the Vulcan CO 2 fossil fuel emission estimates [39,65]. The Vulcan product includes hourly anthropogenic CO 2 emissions estimates at a 0.1 × 0.1 degree (~121 km 2 in the Boston area) spatial resolution for the entire continental US for the year 2002. Vulcan emissions are partitioned into a variety of sectors, including total, residential, commercial, and on-road emissions. For this analysis, total emissions from all sectors were extracted for a 9 pixel area surrounding each tower site. A 9 pixel area was chosen to ensure that the centroids of the three study areas were within 3 km of the tower locations. Annual, seasonal, daily, and diurnal anthropogenic emissions estimates were calculated over the 9 pixel area (1,089 km 2 areas). At all three sites, the emissions and mixing ratios were analyzed by 45° and 180° sections aligned to the cardinal directions. Only the results from Worcester were significantly variable by sector (and only for the east and west 180° sectors) due to the distinct boundaries between vegetation and more intensively developed areas. These emissions represent estimated local emissions, rather than demand-induced emissions elsewhere. For example, the effect of emissions from air conditioning at these tree sites is difficult to assess since most of the electricity was generated outside of our study areas. In contrast, most space heating emissions were captured in our analysis since 82%-90% of all residential homes in the Worcester and Boston areas used natural gas and fuel oil for heating [66].
To estimate CO 2 flux from human respiration, we created a spatially explicit data layer based on a Massachusetts population density map [29] and an estimate of mean per-capita CO 2 flux from respiration (257 g· C· day −1 per person [67]). The population density map was based on block-level data from the 2010 US Census [57] and allocated the population within a given census block to the residential land area [53] of that census block.
At each tower site, we digitized and calculated the canopy cover within a 1km radius using high resolution QuickBird imagery downloaded from Google Maps [68] and the ImageJ image analysis software [69]. The QuickBird images used for this analysis represent land cover from recent, unspecified dates. The vegetative fraction at each site includes only large woody vegetation, excluding lawns and shrubs. Average annual gross primary production and ecosystem respiration at Harvard Forest was 14.0 and 11.5 mg· C· ha −1 · yr −1 , respectively, between 1992 and 2004 [42]. Since long-term flux tower based estimates of gross primary production are not available for urban areas in our region, we scaled gross primary production and ecosystem respiration linearly based on percent canopy cover, using Harvard Forest as our baseline (97% tree canopy cover). In Boston, a robust relationship between canopy cover and biomass was observed [29], supporting the use of canopy cover as a proxy for biomass and leaf area. These first order biogenic flux estimates were calculated for the 1 km 2 area around each site.

Phenology Timing and Time Series
Remotely sensed surface reflectance indices, such as EVI, are commonly used to characterize vegetation properties based on surface greenness and have been widely used to monitor tree phenology and photosynthetic activity [70,71]. Spring and autumn vegetation phenology at each tower site was estimated using a new algorithm that exploits a time series of data from the 30 m resolution Landsat Thematic Mapper satellite sensor [72]. Only landsat pixels with average "summertime" EVI above 0.6 were included in this analysis so that the effect of building shading and drought prone urban lawns-both of which tend to push EVI below 0.6-could be mitigated. While this method reduces biomass estimates in shaded areas, shaded vegetation comprises a very small percentage of total vegetation. We validated this method in Boston using Bing Maps bird's eye pictometry.
The spring and autumn phenological dates roughly correspond to the timing when leaf lengths reach 25% of their seasonal maximum ("green-up") and 90% coloration ("brown-down"), respectively. Annual phenology dates at each pixel were estimated based on the deviation in Landsat observations in 2011 relative to the 1982-2001 long term average phenology at each pixel. Pixels with low seasonal amplitude in EVI were classified as non-forest and removed from the analysis. We then calculated the median of all retrieved phenology dates across a roughly 3 km × 3 km window (100 × 100 pixels) centered on each study site. The seasonal EVI trend was computed using LOESS with a span of 0.5.
Since Boston, Worcester and Harvard Forest are 5, 70, and 99 km away from the coastline, respectively, the Atlantic Ocean likely moderated climate and influenced phenology along our urban to rural gradient. However, given the results of  which highlight clear urban heat islands in cities within 75km of the coast, we feel that we are seeing a primarily urban signal [73].

Spatial Analysis
We examined associations between land cover and CO 2 mixing ratios for 2 different sectors around the Worcester tower site. Each sector spanned 180 degrees, with the "East" sector extending from 0 to 180 degrees, with 0 degrees being polar north. A categorical variable representing the wind sector was then appended to the CO 2 dataset according to observed wind direction, so that each CO 2 observation was associated with one of the two sectors. Since the wind direction observations at Worcester and Boston could have been adversely affected by micrometeorology, they were validated by data from the NOAA affiliated Boston Logan International and Worcester Regional Airport weather stations. For our land cover analysis, we used MassGIS' 2005 impervious surface and land use layers [53], which are based on 0.5 m resolution digital orthoimagery and accessor's parcel information, as well as a 2010-2011 30 m summertime cloud-free EVI layer generated from LandSat Thematic Mapper data [50]. The ISA layer defined an impervious surface as one covered by buildings, parking lots, roads, and compacted, man-made soils. The land cover layer originally contained 40 categories, which were combined into 8 categories including wetland, water, forest, residential, agricultural, open and successional, commercial and industrial, and other developed (including roads). We calculated ISA, land cover, and EVI within a 5km radius of the tower site for each sector (Figure 7) using the ArcGIS 10.0 software [74]. To address a finer spatial domain around the Worcester tower site, we repeated this analysis for 8 sectors representing 45 degree wind fields around the tower (e.g., 0 to 45 degrees, 45 to 90 degrees, etc.).
We also used land cover to determinate a first order estimate of the CO 2 mixing ratio source area at each site. For example, within the 9 × 9 pixel area seen in Figure 2, land cover around the Worcester tower appeared to be dominated by vegetation, despite the presence of downtown Worcester. Yet, there was a statically significant difference in CO 2 mixing ratios between easterly and westerly sectors at the Worcester tower site, presumably due to local emissions sources. Consequently, we assumed that the source area at Worcester (and those at the Harvard Forest and Boston sites) fell primarily within the 5 km radius around each tower since it was this approximate spatial scale that clearly highlighted site specific differences in land cover.

Conclusions
In this study we examined the spatial and temporal variations in atmospheric CO 2 mixing ratios and carbon fluxes across Boston's urbanization gradient. There were large differences in estimated biogenic and anthropogenic carbon fluxes across this gradient with total anthropogenic emissions ranging from 37.3 mg· C· ha −1 · yr −1 in urban Boston to 1.5 mg· C· ha −1 · yr −1 at the rural Harvard Forest. Despite the ~25-fold difference in local emissions, the mean annual difference in atmospheric CO 2 mixing ratios was only 20.6 ± 0.4 ppm from the most rural to the most urban ends of our Boston gradient. The atmospheric signal from vegetation was clear in the observed seasonality at all sites, regardless of the amount of local vegetation (Figure 1). We observed significant differences in growing season length across the gradient, with Boston's growing season exceeding Harvard Forest's by 31 days in 2011. This extended growing season in Boston could potentially increase Boston's annual net carbon sequestration by as much as 50%. In densely populated urban Boston, we estimated that human respiration contributed nearly as much CO 2 as ecosystem respiration with each contributing 2.8 and 3.1 mg· C· ha −1 · yr −1 , respectively, similar to results found in Vancouver, CA, USA [75]. Heterogeneity in land cover across the urban-rural gradient, combined with variable meteorology and concomitant shifts in source areas, created a profound challenge in disaggregating the contributions to the CO 2 signal.
The estimated carbon fluxes in this study highlight potential carbon mitigation strategies in urban areas, which are responsible-directly or indirectly-for the majority of anthropogenic emissions. For example, it has already been suggested that net carbon uptake from vegetation and corresponding soils can be significant, even in urban areas [35]. Woody vegetation in Boston is not actively managed to reduce carbon emissions, but our first order estimate of urban vegetation sequestration is nearly 0.7 mg· C· ha −1 · yr −1 . Grow Boston Greener, which seeks to increase the city's tree canopy cover to 35%, might with time increase the sequestration by an additional 0.2 mg· C· ha −1 · yr −1 (assuming productivity scales linearly with canopy cover), but the establishment and maintenance of those new trees can also have significant associated emissions [36]. Our results suggest other strategies may have a larger effect on the area's carbon budget. For example, urban growth strategies which focus on densifying suburban areas rather than clearing new exurban areas could reduce both transportation related emissions (currently estimated at 10.6 mg· C· ha −1 · yr −1 in Boston) [76] and emissions associated with forest cover loss [30]. Further, residential emissions are the second largest CO 2 source in Boston (8 mg· C· ha −1 · yr −1 ; Figure 2). Given that over 60% of Boston's housing stock was built before 1939 [77], there are ample opportunities for efficiency improvements from this sector [66].
As cities, regions, and nations move forward with climate action plans and treaties, we need to improve our capacity to measure emissions and carbon exchange within urban areas. With our current measurement networks, modeling techniques, and fundamental understanding of the urban carbon cycle, we are not yet able to adequately characterize urban carbon sources and sinks or to define the influence of urban ecosystems on regional atmospheric composition. The significant overlap of biogenic and anthropogenic processes in these areas makes partitioning of atmospheric mixing ratios into component fluxes difficult. Despite scientific uncertainties, current literature clearly suggests pathways for policymakers and highlights the imperative for CO 2 emission reductions. For example, Massachusetts Governor Deval Patrick signed the Global Warming Solutions Act into law in 2008, committing Massachusetts to an aggressive and sustained reduction of greenhouse gas emissions over the next 40 years. To meet these goals, the Commonwealth of Massachusetts and the City of Boston have undertaken a serious effort to develop emissions reduction plans, but current uncertainties in emissions estimates are large and can inhibit effective policy action [78]. To move forward with cost effective and verifiable emission reduction plans, we need to expand our surface observations in urban areas, improve the spatial and temporal resolution of emission estimates, and continue the development of atmospheric models that can integrate such data to provide transparent estimates of spatiotemporal changes in carbon sources and sinks. These advances would provide local policymakers the tools necessary to target emission hotspots and monitor the effectiveness of greenhouse gas emission reduction strategies over space and time.