Urban warming advances spring phenology but reduces the response of phenology to temperature in the conterminous United States

Across the selected 85 cities, the SOS was 6.1 ± 6.3 d (mean ± SD) earlier in urban areas (SOS: 71.3 ± 13.1 d of year) than in corresponding rural areas (SOS: 77.4 ± 15.2 d of year) with a slope of 1.06 (P < 0.001; Fig. 1A) during 2001–2014. In 87% of the study cities (74 cities), the 14-y mean SOS in urban areas was notably earlier than in rural SOS (P < 0.05; SI Appendix, Fig. S2). Significantly negative SOS correlations with preseason temperature (but not with precipitation or shortwave radiation) were identified in all study cities (P < 0.05; Fig. 1B and SI Appendix, Figs. S3 and S4A). Moreover, the $R_T$ was 0.03 ± 0.07 weaker (i.e., less negative) in urban areas, compared with rural ones (slope = 0.85, P < 0.001; Fig. 1B). Of the study cities, 51% (43 cities) further demonstrated significantly weaker urban $R_T$ than the rural ones (P < 0.05; SI Appendix, Fig. S4B). Similar results were obtained after removing the influence of extreme years on $R_T$ (urban $R_T$ was, on average, 0.04 ± 0.15 weaker than rural $R_T$; SI Appendix, Fig. S5) using a temporal bootstrapping approach. The $R_T$ with different preseasons (i.e., January to March, February to March, March, and March to April) was also examined, and the results showed similar weaker urban $R_T$ although with different amplitudes (i.e., 0.02 ± 0.1, 0.02 ± 0.1, 0.03 ± 0.1, and 0.05 ± 0.1, respectively; SI Appendix, Fig. S6).

The SOS significantly advanced with an increase in annual mean air temperature ($T_{avg}$, slope = −1.98 d/°C, P < 0.0001) but not with annual precipitation $\left(P_{rcp}\right)$ (Fig. 2A and SI Appendix, Figs. S7 A and C and S8A and Table S1). However, the magnitude of $R_T$ significantly decreased (less negative) with the advancement of SOS (slope = −0.005/d, P < 0.0001), increase in $T_{avg}$ (slope = 0.009/°C, P < 0.001), and decrease in $P_{rcp}$ (slope = −0.0008/mm, P < 0.01; Fig. 2A and SI Appendix, Figs. S7 A and C and S8 A and C and Table S1). The urban−rural difference in $R_T$ $\left(\mathrm{\Delta }{R}_T\right)$ significantly negatively correlated with the urban−rural difference in SOS (ΔSOS, slope = −0.015/d, P < 0.001; Fig. 2B and SI Appendix, Table S1). In other words, with the increase in the magnitude of urban−rural difference in SOS, the magnitude of urban−rural difference in $R_T$ increased for cities in quadrant II and decreased for cities in quadrant III. Positive Δ$R_T$ (i.e., reduced magnitude in urban $R_T$) was found mainly in cities associated with relatively low $T_{avg}$, and conversely negative Δ$R_T$ (i.e., enhanced magnitude in urban $R_T$) was predominantly found in cities with high $T_{avg}$ (Fig. 2B and SI Appendix, Fig. S8D). The $T_{avg}$ threshold distinguishing the cities with positive Δ$R_T$ from those with negative Δ$R_T$ was 17.3 °C, which captured the sign of Δ$R_T$ in 88% of cities in quadrant II and 72% of cities in quadrant III, respectively (SI Appendix, Fig. S9). This $T_{avg}$ threshold was further used below to differentiate phenological responses to UHI in cold versus warm regions. Geographically, the cities showing advanced SOS but reduced $R_T$ magnitude (quadrant II) were mainly in the northeastern temperate and boreal regions (e.g., Minnesota, Michigan, and Pennsylvania; SI Appendix, Fig. S8 B and D), while the cities showing advanced SOS and enhanced magnitude in $R_T$ (quadrant III) were largely in the southwestern or coastal regions (e.g., Texas, Louisiana, and Florida; SI Appendix, Fig. S8 B and D).

In addition to background climates, the modified local climatic conditions characterized by UHI were closely associated with the magnitudes of change in SOS and $R_T$. All of the cities except two (98%) showed significant UHI effects (ΔT, i.e., T_Urban – T_Rural, 0.6 ± 0.7 °C, P < 0.05; SI Appendix, Fig. S7E). However, the responses of ΔSOS and $\mathrm{\Delta }{R}_T$ to ΔT differed from cold regions to warm regions separated by 17.3 °C as the threshold (see above). For cities in cold regions, ΔSOS and Δ$R_T$ significantly decreased at −7.0 d/°C (P < 0.01; Fig. 3A) and increased at 0.04/°C (P = 0.05; Fig. 3C), respectively, with the intensification of UHI. In warm regions, the intensification of UHI significantly strengthened the urban $R_T$ reduction at 0.07/°C (P < 0.05; Fig. 3D) but not SOS (P = 0.06; Fig. 3B). It’s worth noting that Δ$R_T$ was negative when ΔT was less than 0.5 °C (Fig. 3D). Such enhanced urban SOS response to temperature gradually changed to a reduced one (positive Δ$R_T$) with the intensification of UHI in warm regions. These contrasting patterns show that the ΔSOS and Δ$R_T$ widened in cold regions, but more complex changes occurred in warm regions with UHI intensification.

Two process-oriented phenology models, namely, the alternating (ALT) (26) and parallel (PAR) (27) models, were used to examine the possible causes underlying the SOS and $R_T$ changes in cities. These two models, which have been shown to agree well with both in situ and remote sensing observations, require daily air temperature as the only input to predict spring budburst (28). They assume that plants accumulate chilling and forcing during winter and spring, and spring budburst occurs when forcing and chilling accumulations meet certain requirements (SI Appendix). Both models captured the interannual variation and spatial variability of SOS with a correlation of 0.85 between observation and prediction during the independent validation for 2012–2014 (P < 0.001; SI Appendix, Table S2 and Fig. S10). The models also reproduced the ΔSOS as observed, albeit with underestimated magnitudes (SI Appendix, Fig. S11). Compared with those in rural areas, urban plants were simulated to accumulate more forcing and less chilling for 78% (ALT) and 67% (PAR) of cities and years (quadrant II; SI Appendix, Fig. S12). Similar results were obtained when comparing 14-y mean forcing and chilling accumulations in urban versus rural areas in 85 cities; that is, 85% (ALT) and 67% (PAR) of cities showed more forcing and less chilling accumulations in urban areas (quadrant II; SI Appendix, Fig. S13). How the changes in forcing and chilling account for the changes in SOS and $R_T$ is discussed below.

This study provides direct evidence that urban warming not only advanced SOS but reduced the magnitude of SOS response to temperature in 85 large US cities. Such urbanization-induced phenological changes are broadly consistent with recently reported research on natural ecosystems (21, 22, 29). As a preview of possible impacts of future warming on phenology, this finding suggests that, while spring phenology will likely continue to advance in natural ecosystems, the rate of SOS advancement will slow down under anticipated climate warming conditions (6, 7, 21, 22). To explain the reduced $R_T$ magnitudes, we consider four hypotheses: 1) thermal budget percentage (6), 2) insufficient chilling (7, 30), 3) reduced preseason length (31), and 4) photoperiod restriction (32).

The hypothesis of thermal budget percentage suggests that spring phenology is more sensitive to temperature in colder, higher-latitude sites than that in warmer regions because small absolute changes in temperature constitute greater relative changes in thermal balance under colder conditions. This theory thus justifies the reduced magnitudes of urban $R_T$. If we treat urban and rural areas as warm and cold conditions, respectively, rural SOS is more responsive to warming because a smaller increase in temperature constitutes a larger thermal balance in rural areas, leading to a relatively larger magnitude for rural $R_T$ than that for urban $R_T$ (i.e., positive Δ$R_T$).

The hypothesis of insufficient chilling assumes that, if the necessary chilling for dormancy break is not fully met (e.g., under warm winters), plants become less responsive to spring warming. Specifically, after growth cessation in early fall, plants need to be exposed to low chilling temperatures for a certain period to break dormancy. After that, phenological development (e.g., budburst) occurs as a result of high forcing temperatures above a threshold (33). If the necessary chilling is not fully met, the forcing process might be delayed or the forcing requirement might be increased, leading to a decreased magnitude in the temperature response of SOS. According to this hypothesis, the chilling requirement may not be fully achieved in urban areas, due to mild winters, causing an increased requirement for forcing and, taken together, reduced magnitude in temperature response of SOS in urban areas (i.e., less negative $R_T$ or positive Δ$R_T$). Simulations from the two phenology models also support this hypothesis (SI Appendix, Figs. S12 and S13).

The hypothesis of preseason length suggests climate warming reduces the magnitude of $R_T$ simply by reducing the length of the preseason, due to a faster progression toward budburst in spring (31). However, we did not find significant changes (P = 0.18) in preseason length from urban (2.0 ± 0.8 mo) to rural areas (1.9 ± 0.7 mo; SI Appendix, Fig. S14). Thus, the impacts of preseason length on $R_T$ can be excluded in this study.

Photoperiod (i.e., the period that a plant is exposed to light in a day) has been proposed to restrict the advancement of phenology under warming by imposing a hard limit (5, 22, 34). However, we did not find evidence that a hard limit had been reached in this study. Photoperiod is a cue for the progression of phenology at a given location because it gradually increases from winter to summer and then decreases until the following winter (35). Unlike large fluctuations in temperature, photoperiod does not vary from year to year, and thus has a stabilizing effect on phenology (36). Experimental and modeling studies revealed that, although a warmer spring leads to an earlier budburst, the resulting shorter day length when budburst occurs might have an adverse effect on phenology by slowing down cell growth in the spring to protect leaves against frost damage (32). In this study, urban and rural areas were in close proximity and experienced nearly identical photoperiods on any given day under natural conditions; however, the urban SOS actually occurred earlier at times with progressively shorter day lengths. Photoperiod could therefore partially contribute to the reduced $R_T$ magnitudes along with the insufficient chilling restriction. However, it is likely that the hard limit of photoperiod was not reached in this study, since SOS still showed a continuous advancement in warmer conditions.

Most cities in the cold region showed reduced $R_T$ magnitudes, but many in the warm region showed enhanced $R_T$ magnitudes (negative Δ$R_T$; quadrant III in Fig. 2B). Cities associated with enhanced $R_T$ magnitudes were mainly located in southern or western coastal regions (e.g., Texas, Louisiana, and Florida; SI Appendix, Fig. S8D), where winters are usually mild and short, followed by hot and long summers. These cities also exhibited complex phenological responses to UHI effects compared with those in the cold regions (Fig. 3 B and D). Such complex phenological responses are possibly caused by the dependence of local vegetation on other environmental cues for spring green-up in addition to temperature. For example, Prunus dulcis (sweet almond), a species native to Iran but widely cultivated elsewhere (e.g., California and Utah), is particularly sensitive to the available soil water (37).

Both phenology models are representative of the forcing and chilling processes and are consistent with the second hypothesis discussed above. They captured most of the observed SOS variations over time and among cities (SI Appendix, Fig. S10) and, more importantly, the observed earlier urban SOS (SI Appendix, Fig. S11). As these models are dependent on temperature only, and although other factors may influence the urban phenology to a certain extent (discussed below), urban warming dominated the phenological changes in this study. Because forcing processes determine how fast buds develop in the spring (5), a more and faster accumulated forcing, as estimated in the urban areas (SI Appendix, Fig. S12), explains the earlier urban SOS (Fig. 1A). In contrast, the increased accumulated chilling reduces the forcing requirement for budburst (5), as described in both models (SI Appendix). Thus, with less chilling accumulated (SI Appendix, Fig. S12), urban plants require more forcing for budburst, which offsets the forcing-induced early SOS, to a certain extent, and accounts for the decreased $R_T$ magnitudes (38) (Fig. 1B).

The modeled Δ forcing and Δ chilling (SI Appendix, Figs. S12 and S13) also explained the widening gaps of urban−rural difference in SOS and $R_T$ with intensified UHI for cold regions (Fig. 3 A and C). This is because the intensified UHI can cause continuous forcing gain and chilling loss for urban plants and thus further advance SOS and reduce $R_T$ magnitudes. For warm regions, however, the ΔSOS and Δ$R_T$ pattern (Fig. 3 B and D) could not be explained by such temperature-dependent phenology models (SI Appendix, Fig. S12). This is likely due to the complex response of spring phenology to various environmental cues as discussed above. In addition, the slopes of Δ forcing against Δ chilling in warm conditions were steeper than those in cold conditions (SI Appendix, Figs. S12 and S13). This suggests that, with the same amount of chilling loss in urban areas, the gained forcing is greater in a warm region than that in a cold region. Such differences in warm versus cold regions were mainly due to the differences in the relative proportion of forcing and chilling accumulations. For example, in cold regions, the chilling period is longer, and, consequently, the chilling accumulation is greater than that in warm regions. Therefore, with the same amount of chilling loss in urban areas, the relative portion of chilling lost is smaller in cold regions, causing a smaller forcing gain.

The influence of urban size (an indicator of urbanization level) on vegetation phenology was investigated, but no significant relationship was found between the urban size and ΔSOS or Δ$R_T$ (SI Appendix, Table S1). A possible explanation is that various background climate conditions and associated plant functional types could cause different phenological responses to urban warming, even for cities of similar size (39). Other abiotic and biotic factors that were not considered in this analysis (e.g., urbanization-induced changes in nitrogen deposition, precipitation regimes, and plant species) could also affect urban phenology differently from temperature (9). For example, spring green-up dates vary greatly among different plant species. Because urbanization has partly altered the species and composition of vegetation in urban ecosystems (e.g., invasive species Acer platanoides in northeastern cities and the tropical species Jacaranda in California), diverse plant species may cause highly spatially heterogeneous phenology in cities (40). We nonetheless believe that our broader argument is valid. First, despite exotic species in cities, previous studies found that the rich variety of tree species, native or exotic, does not change significantly with urbanization gradient (41). With a considerable amount of native species in cities, such mixed phenological signals still reflect the phenology of tree communities in their geographic regions (e.g., rural areas) (42). Second, exotic species in cities have diverse phenology niches, some of which are earlier and others are later than those of native species (43). Because the employed satellite observations were derived from a mixed signal of all species within a grid at 500 m by 500 m spatial resolution, the exotic species could not account for the consistent SOS changes across most cities as documented in this study.

In summary, we reported an urban warming-induced advanced SOS and weakened SOS response to temperature change in 85 cities across the United States and demonstrated that increased forcing and reduced chilling were largely responsible for these phenological changes. These findings provide observational evidence of a weakened temperature response of SOS (i.e., decreased $R_T$ magnitude) under urban warming conditions that deserves attention, along with the widely acknowledged earlier urban SOS. Moreover, such reductions in $R_T$ magnitude were found to be mainly associated with cities in cold regions, and their changes correlated well with the intensity of UHI, whereas disparate $R_T$ patterns existed in warm regions. These results suggest that the previously identified advancement in spring phenology in natural ecosystems may continue but that the rate of advancement will likely slow down under future warming, especially in cold areas.

The spring phenology indicator SOS during 2001–2014 was obtained from a MODIS Land Cover Dynamics (MCD12Q2) product (23). The climate data during 2000–2014 were obtained from the TopoWx dataset (24) and the Daymet dataset (25). The National Land Cover Database was used to classify pixels of urban areas, deciduous forest, evergreen forest, and mixed forest in this study (44). The urban extent was derived from the Defense Meteorological Satellite Program/Operational Linescan System nighttime light data (45). Cities with urban areas larger than 500 km² were selected in this study (46).

To understand the main driving climate factor(s) of SOS, covariations (i.e., partial correlation coefficients) between SOS and climate variables (i.e., $T_{avg}$, precipitation, and shortwave radiation) during 2001–2014 were calculated using partial correlation analysis. To estimate the uncertainty of the $R_T$ obtained from the partial correlation analysis, we employed a nonparametric bootstrapping method (SI Appendix, Materials).

This work was supported by the Terrestrial Ecosystem Science Scientific Focus Area (TES SFA) project, the Reducing Uncertainties in Biogeochemical Interactions through Synthesis and Computing Scientific Focus Area (RUBISCO SFA) project, and the Energy Exascale Earth System Model (E3SM) project, funded by the US Department of Energy (DOE), Office of Science, Office of Biological and Environmental Research. Oak Ridge National Laboratory is supported by the Office of Science of the DOE under Contract DE-AC05-00OR22725. L.M. is also supported by NASA Research Opportunities in Space and Earth Science (ROSES) Climate Indicators and Data Products for Future National Climate Assessments (INCA) Program “NNH14ZDA001N-INCA” and Future Investigators in NASA Earth and Space Science and Technology (FINESST) Program “NNH19ZDA005K.” A.D.R. acknowledges funding from NSF, through the Macrosystems Biology program (DEB-1702697). This manuscript has been authored by UT-Battelle LLC under Contract DE-AC05-00OR22725 with the DOE. The US Government retains and the publisher, by accepting the article for publication, acknowledges that the US Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US Government purposes.