Disentangling the trend in the warming of urban areas into global and local factors.

Large cities account for a significant share of national population and wealth, and exert high pressure on local and regional resources, exacerbating socioenvironmental risks. The replacement of natural landscapes with higher heat capacity materials because of urbanization and anthropogenic waste heat are some of the factors contributing to local climate change caused by the urban heat island (UHI) effect. Because of synergistic effects, local climate change can exacerbate the impacts of global warming in cities. Disentangling the contributions to warming in cities from global and local drivers can help to understand their relative importance and guide local adaptation policies. The canopy UHI intensity is commonly approximated by the difference between temperatures within cities and the surrounding areas. We present a complementary approach that applies the concept of common trends to extract the global contributions to observed warming in cities and to obtain a residual warming trend caused by local and regional factors. Once the effects of global drivers are removed, common features appear in cities' temperatures in the eastern part of the United States. Most cities experienced higher warming than that attributable to global climate change, and some shared a period of rapid warming during urban sprawl in the mid-20th century in the United States.

Introduction

Three major anthropogenic factors and their interactions influence climate at the city scale: the effects on global climate of the accumulation of radiative active substances in the atmosphere; the impacts on regional climate of air pollutants, such as aerosols and black carbon; and the urban heat island (UHI) produced by urbanization, which entails the replacement of natural landscapes by denser, higher thermal capacity materials and structures, the generation of urban canyons, and higher anthropogenic heat release. , , Urban climate is also determined by natural factors. Some are static (i.e., time invariant), such as geographical location and altitude, while others are dynamic, such as the Atlantic Multidecadal Oscillation (AMO) and other regional and large‐scale climate variability processes.

The UHI effect is given by the difference in the warming of a city with respect to its preurban conditions. It can be explained in terms of surface energy balance processes related to short‐ and long‐wave radiation exchange, and latent, sensitive, and conductive heat flows. , The magnitude of the UHI is a function of population, urban morphology, and physical characteristics, such as city size, compactness, gross building volume, and anisometry (city shape), as well as vegetation and waste heat release. , , , , , In the case of large cities, UHI warming can be similar in magnitude (about 2.5–4.5 °C) to that expected under high emissions climate change scenarios for the end of this century. , , The UHI is commonly characterized by great spatial variability, large differences in structure and intensity, and possibly multiple UHIs across a single city. UHI intensities tend to be higher near the downtown area and to decrease toward the suburbs and urban periphery. The UHI varies greatly between seasons and throughout the day, being typically stronger at night. The three types of UHI commonly considered in the literature are the surface, canopy, and boundary layer UHI. Remote sensing is typically used to investigate the surface UHI, while measurement of the boundary UHI is based on sensors installed in towers, tethered balloons, or radiosondes. The canopy UHI is characterized by comparing air temperature records from weather stations located within a city and its rural surroundings. , , Some approaches to approximate the UHI intensity include empirical functions relating population counts and urban warming, , , using urbanization indicators, comparing the values of the regional and local trend slope estimates, contrasting surface temperatures from reanalysis reconstructions , (e.g., climate model simulations that do not account for land‐use differences so that the differences between observations and model simulations should be caused by the UHI), analyzing climate model simulations, , and using satellite data to derive urban land surface temperatures. , The proposed method and the analysis presented here refer to the canopy UHI.

Several studies found that the logarithm of a city's population is positively correlated with the UHI intensity such that more populous cities experience larger changes in local climate. , , This relationship can be modulated by the background climate and other factors, in particular by the total amount of precipitation. Urban sprawl is a typical mode of expansion and can lead to a rapid expansion of the UHI over the city, although characteristics, such as urban design, with low compactness ratio can mitigate it. , By contrast, compact cities may lead to a more intense UHI that extends over a smaller area. Urban sprawl has also been associated with increases in annual extreme heat events.

The U.S. population increased more than 40‐fold since 1800, from about 5 million to over 300 million in 2010, implying a significant transformation of the physical landscape and ecological systems. , About 80% of the U.S. population lived in urban areas in 2010 that covered only 3.1% of the total land area. While urban development started in the northern and coastal southern states and then spread to other parts of the country, the peak of the expansion and densification in most regions was reached in the period 1950–1975. Population is highly correlated with urban built area, and population growth has followed the urbanization expansion and densification patterns. These changes were accompanied since the 1920s by urban sprawl, which reached its period of most rapid increase after WWII and continued until the mid‐1990s.

The literature on detection and attribution of urban warming is rapidly expanding, in part motivated by the potential bias that the UHI effect in cities can impart to average global and regional temperatures. Also of interest is the goal to quantify the UHI effect and to inform climate adaptation policy at the urban scale. , Efforts have been devoted to evaluate the magnitude of the bias the UHI could impart to the global warming trend and to devise ways to remove it. , The effect of the UHI on local temperatures is undisputed, but there is no consensus about the extent of its influence on regional and large‐scale temperatures. , The contribution of the Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) concluded that the UHI effect on temperatures is negligible on large‐scale trends. For global temperature, it is unlikely that it could represent more than 10% of the centennial warming trend. However, the influence of the UHI could be much larger on regional temperatures. Recently, the study of the effects of the UHI on regional temperatures has focused on China, given its rapid development and high rates of urbanization. Some studies have found that the effect of UHI on the warming trend at the national level could be up to 0.14 °C per decade, although studies using other methods suggest a more modest influence on the warming trend of 0.03–0.05 °C per decade. For Northern China, an urban warming trend of 0.11 °C per decade was detected, which accounts for about 37.9% of the total warming in the region. , While the UHI has a large impact on local temperatures, in the case of area‐weighted annual mean temperatures for the continental United States, its effect has been deemed by some authors as negligible (about 0.06 °C), , although others have found this effect to be considerable, contributing up to 0.27 °C per century. At the local and city scales, the UHI effects are much larger and can account for significant shares of the total warming trend. In New York City, about a third of the warming experienced since 1900 has been attributed to the UHI intensification. This warming is spatially heterogeneous, and not easily explained in relation to land surface characteristics, as exemplified by the similarity between the rate of warming in Central Park and that in other parts of the city. In the case of cities in developing countries, the share of warming explained by the UHI effect can be much larger, for example, in Beijing and Wuhan, where urban warming accounts for about 65–80% of the total warming over the 1961–2000 period.

In contrast with some previous studies, we focus on removing the effect (bias) in local temperatures that can be attributed to global climate change in order to produce a residual trend that is mainly imparted by local/regional forcing factors. Hence, our approach is similar in spirit to a signal extraction problem. For this purpose, we draw on recently proposed time‐series methods for the attribution of global climate change based on the existence of common trends between globally aggregated radiative forcing and temperature series. The attribution of global climate change has been an active topic of research during the last decades and has shown the existence of a clear signal that can be traced back to globally aggregated anthropogenic forcing in a variety of climate variables, from global to regional scales. , , , , However, this type of attribution studies at scales finer than hemispheric are based almost exclusively on the results from physical climate models, with very few exceptions. Here, we focus on extending current observations‐based global climate change attribution methods to purge the global climate change signal from urban annual surface air temperatures and to characterize in more depth the residual warming trend. We use UHI detection and attribution methods to investigate the possible contribution of local forcing factors to the residual warming trend. This is relevant from a climate research perspective as it allows a better understanding of the contributions of urbanization processes and of the effects of increases in global radiative forcing factors at the local scales. It is also relevant from a decision‐making perspective as it can help to define city and national mitigation and adaptation policies. Recent studies have stressed the importance of assessing local and global climate change in defining the economic impacts in cities, in which more than half of the global population lives and about of 80% of global GDP is produced.

Methods and data

Attribution of the global warming trend at the city level and estimation of the residual trend due to local and regional forcing factors

Commonly used procedures tailored to estimate the magnitude of the warming caused by the canopy UHI effect are based on the differences between temperatures for stations located within and outside the city. However, this does not allow to purge the global warming signal from that caused by local and regional factors that contribute to the observed warming, nor to extract features of its temporal evolution. Moreover, it does not consider the effects of large‐scale natural variability. The selection of rural stations was shown to be complex and potentially has nontrivial effects on the resulting estimates. Also, the comparison between rural and urban stations can introduce an unnecessary noise component present in temperatures from weather stations outside of the city. Our proposed methodology addresses some of these problems and is easy to implement. It also extends the literature on the attribution of global climate change to the urban/local scales and connects it to local‐scale processes, such as the UHI.

We focus on two characteristics of the warming trend to separate the effects of local and global factors in city‐level temperatures: the magnitude of the warming and the shape of the warming trend. These characteristics are investigated by comparing local‐scale (city) and large‐scale temperatures series obtained from wider regions containing the selected cities.

We first investigate the contribution of global factors to the warming trend of city‐level temperatures by testing for the existence of a common trend via a cotrending test applied to regional and local temperatures and some measures of aggregate total and anthropogenic radiative forcing series. The existence of a common trend suggests that at least part of the secular movement contained in regional and local temperatures is produced by global‐scale changes in the Earth's energy balance. Hence, part of the observed warming at the different scales is then attributable to globally aggregated anthropogenic forcing factors.

For the local temperature series, this common trend can be modulated by the superimposition of other warming trends caused by the local and regional effects of external forcing factors, including urbanization and land‐use change. Observed regional and local temperatures can differ significantly over some periods of time and this can be explained by the evolution of relevant local and regional determinants.

We shall present a regression‐based method to isolate and analyze the features of the warming trends in local temperature series. To illustrate the proposed method, we analyze 30 cities and discuss the similarities and differences about how the evolution of local factors, such as urbanization, could have affected urban climate within some geographical regions. Consider the following calculation of regional average temperature: where are area weights and is the area‐weighted regional temperature and (i = 1,…,n) the local temperature for location (or grid cell) i. If for all i (all grid cells have the same area) and ω is normalized to 1, then (1) simplifies to . As with any time series, can also be decomposed to trend (), cycle (), and irregular () components:

Focusing first on the trend component, we have: where n is the total number of temperature series in the selected region, k and m are the number of urban and nonurban locations, respectively, and . The selected region should be large enough that n is much larger than k. Note that while all temperature series in the selected region contain a global warming trend, city temperatures can have both local and global warming trends and thus the total number of trends is . This decomposition can be seen as a modified version of Lowry's additive scheme in which both urban and rural temperatures are a function of background climate (represented here by the common global warming trend), the effects of local factors, and, for cities, the effects of urbanization. The trend is a weighted average of the warming trends due to urban (local) factors and to large‐scale climate change , which is imparted by changes in the globally averaged total radiative forcing (TRF). Since the global warming trend is common to all temperature series, the second term in the right‐hand side of Eq. (3) becomes . In the case of the local warming trends, there are at most common trends (note that the normalizing factors are such that ). However, for this condition to occur, all urban locations should share the same warming trend. This is unlikely since the drivers of the local warming trend (i.e., urbanization, land‐use change, and aerosols emissions, among others) are quite heterogeneous so that the local trends are expected to have unique features. Hence, the temperature series for each city contains the global warming trend , which is identical for all cities since it is imparted by globally aggregated radiative forcing. However, in addition, cities experience a local warming trend caused by UHI and other local factors, denoted by , which is city‐specific, as it responds to a particular combination of local drivers (e.g., population, city shape, vegetation, water bodies, and pollution).

With local warming trends being heterogeneous and each weighted by a normalizing factor equal to , the local trends get diluted, while the global warming trend becomes clearer at the regional scale due to averaging. As m gets larger than k, urban warming has less influence on the magnitude of large‐scale warming and will not alter the shape of the global warming trend present in large‐scale temperatures, nor modify its features. In this case, . As mentioned in the introduction, in practice, some residual effect of UHI can still be present and this representation is only approximate, though as we shall document still very useful.

The effects of natural variability oscillations of large‐scale processes, that is, the cyclical components , are expected to be common to all local temperature series within a selected region, while idiosyncratic local variations, or local cyclical components , are unique or shared by a limited number of locations. Averaging produces a clearer estimation of the common cycles and reduces the effects of idiosyncratic local variations. Consider the aggregate cyclical component:

For simplicity, we assume that all n series share a common cycle and d of them have an additional unique cyclical component. Then, the common cycle is weighted by the normalizing factor and the idiosyncratic cycles by ; in the case of common local cycles, they would be weighted at most by . As n gets larger than d, the cyclical component of is increasingly dominated by the common large‐scale variability component; that is, and the effects of local variability are negligible. As in the previous case, this representation is approximate as some residual regional oscillations can be present, though still useful.

Owing to the existence of a common warming trend, city temperatures () can be represented as a function of regional temperatures plus local trend and cyclical components :

If urban and regional temperatures warm at the same rate over long periods of time, then . Imposing this coefficient restriction and subtracting from both sides of Eq. (5), we obtain:

Using Eq. (6), the common trend and cycle produced by large‐scale climate variability cancel out and the idiosyncratic trend and cyclical components can be estimated. The effects of differences in latitude and elevation could affect the assumption that . We briefly discuss this issue. The observed response of annual temperatures to changes in external forcing was recently studied, using similar methods as those used here, showing that the transient response is only statistically different from the other regions for the high latitudes in the Northern Hemisphere, because of the Arctic amplification phenomenon. , While absolute temperature values can show large spatial variability, the change in temperatures due to modifications in globally aggregated radiative forcing tends to be more spatially homogeneous. , Some studies suggest that elevation may amplify the rate of warming, although these results are mixed and some other studies have found no clear correlations, arguing that the association may be due to inhomogeneities in the data. , ,

Once is obtained, the magnitude of the local warming component can be estimated, its features can be described (i.e., identification of periods of faster/slower warming), and its drivers can be investigated. However, it is important to consider that in practice, may be contaminated by the influence of regional factors not accounted for by the large‐scale warming trend and thus can be an imperfect representation of the warming trend generated by local factors. As such, the residual warming trend can be expressed as , where is the true, unobserved warming trend due to the UHI effect, and is a disturbance component that can include the effect of multiple local factors other than UHI (e.g., elevation and vegetation), as well as regional trends that are not included in and a potential UHI bias in large‐scale temperature. As mentioned in the introduction, when large‐scale temperatures are considered, this effect is likely small but it can depend on the region. ,

In cities for which no local processes contribute to the urban warming trend, would be a stationary process, as the common trend would cancel. We can test for remaining trends or other features in that would indicate additional effects of local factors on local temperatures. We focus on two complementary representations of the nonstationarities in that provide useful information about the effects of local factors, which can help to identify their drivers. The first is the estimate of the slope coefficient B in the following regression: which provides an estimate of how local factors have contributed to making the warming rate at the urban scale higher or lower across different cities compared with the large‐scale warming trend. Second, can include nonlinearities, such as level shifts. Estimates of the level shifts and their break dates can be obtained using the following linear regression model and applying multiple break tests as described below to estimate the value of the parameters . The results can suggest how local factors have modified the warming trend at the city scale, producing deviations from the common trend. If there are level shifts in , their origin can be investigated using proxy variables (such as population counts) for factors related to local climate, such as urbanization and urban sprawl. Moreover, the dates of the level shifts can be related to the occurrence of other features in variables that influence local climate, such as urbanization and population. The tests for cotrending and structural change are described in the following paragraphs.

Bierens nonparametric nonlinear cotrending test

For the Bierens cotrending test, the common nonlinear trend included in the time series does not have to be parameterized. The vector of series is modeled via the following trend noise decomposition , where , with a k‐variate time series, a k‐variate zero‐mean stationary process, and a deterministic k‐variate general nonlinear trend function that allows for structural changes. If there exists a nonzero vector θ such that , the system of variables is nonlinear cotrending. The null hypothesis of this test is that the multivariate time series is nonlinear cotrended, and one or more linear combinations of the time series are stationary around a constant or a linear trend. Note that this test is a cointegration test when it is applied to series that contain unit roots. The test is based on the generalized eigenvalues of the matrices M 1 and M 2 defined by , where if and if , and are the estimates of the vectors of intercepts and slope parameters in a regression of on a constant and a time trend. Also, , where , with T the number of observations and . Solving and denoting the solution of the largest rth eigenvalue by , the test statistic is . The null hypothesis is that there are r cotrending vectors against the alternative of r–1 cotrending vectors. The test has a nonstandard distribution and the critical values are tabulated in Bierens. The existence of r cotrending vectors in r+1 series indicates the presence of r linear combinations of the series that are stationary around a linear trend so that these series share a single common nonlinear deterministic trend. Such a result indicates a strong secular comovement in the r+1 series.

Bai–Perron multiple structural change test

Bai and Perron , , provided a framework to test for multiple structural changes in which some of the parameters in a regression model are allowed to change at m unknown break dates. Consider the linear regression model with m breaks (m+1 regimes) , where , and are vectors of covariates, β is a vector of fixed coefficients, while is a vector coefficients subject to change across regimes, and is the error component. The number of breaks and the break dates are unknown and estimated by minimizing the sum of squared residuals. Three different types of procedures are used to estimate the breakpoints: (1) global maximizers of breakpoints; (2) sequentially determined breakpoints; and (3) a mixture of global and sequential procedures. Here, we use the sequential L breaks versus L+1 breaks procedure. The critical values for these tests were tabulated by Bai and Perron.

Data

City‐level annual surface air temperature data were obtained from the National Oceanic and Atmospheric Administration (NOAA; NOAA National Centers for Environmental Information, Climate at a Glance: City Time Series, published October 2020, retrieved on October 10, 2020 from https://www.ncdc.noaa.gov/cag/). Although these time series were chosen by NOAA to represent temperatures in cities, a significant portion are in airports which are rarely close to the most urbanized part of the city where the main UHI is expected to be. This is a frequent limitation in studies trying to characterize the UHI effect as the intensity of the UHI is likely underestimated. The region of study extends over the coordinates −90°E to −65°E and 35°N to 50°N, and the average large‐scale temperatures for this domain were obtained from the NOAA GlobalTemp_v5.0 merged land–ocean surface temperature dataset and processed in the KNMI Climate Explorer data portal (https://climexp.knmi.nl/). Note that this region covers an area of about 3.5 million km2 and, as discussed previously, the effect of any single city temperature series on the area‐weighted average is expected to be small or negligible. Temperature records for cities located within this domain were considered if they had no more than two consecutive missing observations. The missing data were filled in using linear interpolation between adjacent observations. A total of 30 cities were then selected for the analysis (Table 1; see Table S1, online only, for station location, name, and details). These cities are located in 18 states: New York, Massachusetts, Maine, Michigan, Pennsylvania, Virginia, Wisconsin, West Virginia, Kentucky, Washington DC, Indiana, Maryland, Illinois, North Carolina, Ohio, Vermont, Tennessee, and New Hampshire. Figure S1 (online only) shows the region of study and the location of each station, while Figures S3–S32 (online only) provide a satellite image of each city and the station position. Most are located in airports and parks. However, some have been integrated into the city or are now located in significantly urbanized areas (e.g., Figs. S1, S2, S13, and S23, online only), while some are in moderately urbanized or almost rural areas (e.g., Figs. S10, S11, and S19, online only). Seven of the chosen cities have observations for the period 1895–2019, while the temperature series for the others have shorter samples (Table 1; Fig. S2, online only). Data on the elevation of the weather station are from https://www.ncdc.noaa.gov/cag/city/data‐info (Table S1, online only). Radiative forcing variables were obtained from the Goddard Institute for Space Studies of NASA (https://data.giss.nasa.gov/modelforce/). The aggregated radiative forcing variables used are the well‐mixed greenhouse gases (WMGHGs: carbon dioxide (CO2), methane (NH4), nitrous oxide (N2O), and chlorofluorocarbons), and the TRF, which includes WMGHGs plus ozone (O3), stratospheric water vapor, solar irradiance, land‐use change, snow albedo, black carbon, reflective tropospheric aerosols, and the indirect effect of aerosols. Figure S2 (online only) shows the regional and city average annual temperatures, as well as TRF and WMGHGs. AMO was obtained from NOAA's Physical Science Laboratory (https://psl.noaa.gov/data/timeseries/AMO/). City population counts (Fig. 1) were obtained from the Spatial History Project of the Center for Spatial and Textual Analysis at Stanford University.

Table 1

Estimates of local warming rates, shifts in mean, and peak populations

	Sample	B  °C/100 years	δ1	δ2	Pop peak
New Yorku (New York)	1895−2019 (125)	1.04**			8,475,980
Bostonu (Massachusetts)	1936−2019 (84)	0.66**	0.36 [1949]	0.27 [1990]	801,444 [1950]
Portlandu (Maine)	1895−2019 (125)	0.63**	0.47 [1947]		77,634 [1950]
Muskegonu (Michigan)	1915−2019 (105)	0.48**	0.43 [1953]		48,429 [1950]
Erieu (Pennsylvania)	1926−2019 (94)	0.74**	0.46 [1959]		138,440 [1960]
Norfolku (Virginia)	1911−2019 (109)	0.21*	0.28 [1989]		307,951 [1970]
Green Bayr (Wisconsin)	1895−2019 (125)	0.25*	0.38 [1987]		104,879
Elkinsr (West Virginia)	1926−2019 (94)	−0.35*			9121 [1950]
Lexingtonr (Kentucky)	1895−2019 (125)	−0.06			323,780
Allentownu (Pennsylvania)	1925−2019 (95)	1.35**	0.45 [1964]	0.47 [1990]	122,623
Buffalou (New York)	1895−2019 (125)	0.05			256,052 [1950]
Albanyu (New York)	1895−2019 (125)	0.33**	0.29 [1949]		96,853 [1950]
Washington Reaganu (Washington DC)	1946−2019 (74)	0.81**	0.23 [1964]	0.29 [1989]	720,687 [1950]
Philadelphiau (Pennsylvania)	1948−2019 (72)	1.08**	0.45 [1990]		2,071,605 [1950]
South Bendu (Indiana)	1948−2019 (72)	0.14			132,445 [1960]
Baltimoreu (Maryland)	1940−2019 (80)	0.53**	0.27 [1985]		949,708 [1950]
Springfieldr (Illinois)	1901−2019 (119)	−0.39**	−0.31 [1974]		116,250
Syracuser (New York)	1939−2019 (81)	0.25			220,583 [1950]
Indianapolisu (Indiana)	1948−2019 (72)	0.36			820,445
Charlotter (North Carolina)	1940−2019 (80)	0.45*			731,424
Milwaukeeu (Wisconsin)	1939−2019 (81)	1.71**	0.86 [1984]		741,324 [1960]
Clevelandu (Ohio)	1939−2019 (81)	0.74**	0.41 [1986]		914,808 [1950]
Sault Ste. Marieu (Michigan)	1931−2019 (89)	0.51*	0.38 [1960]		18,722 [1960]
Burlingtonu (Vermont)	1895−2019 (125)	0.51**	0.49 [1988]		42,417
Molineu (Illinois)	1944−2019 (76)	0.10			46,407 [1980]
Fort Wayner (Indiana)	1941−2019 (79)	0.42*			253,691
Greensborou (North Carolina)	1929−2019 (91)	0.16			269,666
Chattanoogau (Tennessee)	1939−2019 (81)	0.50**			169,514 [1980]
Johnson Cityr (Tennessee)	1938−2019 (82)	−0.48**	−0.32 [1976]		63,152
Concordu (New Hampshire)	1940−2019 (80)	0.52**	0.48 [1957]		42,695

Note: For the “Sample” column, the sample size is given in parentheses. *,** denote statistical significance at the 5% and 10% levels, respectively. Standard errors were calculated using the Newey–West correction (using the default values in Stata). Superscripts r and u denote stations for which the surroundings are mostly rural or urban, respectively. The number of breaks is determined using the Bai and Perron methodology. Columns 4 and 5 present the estimate of the parameters from regression (8). Empty cells under δ1 indicate no breaks, while empty cells under δ2 indicate no second break. The last column shows the population peaks, and the year in which this peak was achieved is shown in brackets.

Figure 1

Population counts per city for the period 1895–2010.

Results and discussion

Cotrending tests were first applied to the set of variables WMGHGs, TRF, and . The results strongly suggest the existence of a single common nonlinear trend between these series. Moreover, by construction, it can be concluded that the trend is imparted from WMGHGs to TRF and to (Table S2, online only, first panel). Linking the cotrending test result to an energy balance model , allows to attribute a dominant part of the observed warming in to anthropogenic forcing. We removed the trend imparted by TRF to approximate the cyclical component of . This was done via the regression and defining . Figure S33 (online only) shows that the low‐frequency cyclical component in closely follows that in AMO.

Visual inspection of Figure S2 (online only) suggests the existence of a possible common trend between radiative forcing and regional and local temperatures, although this trend is masked at the local level by large variance and the effects of local factors that can distort it. Applying the cotrending test to groups of four time series containing WMGHGs, TRF, , and for each city allows to conclude that the same warming trend is also present in temperatures at the city level (Table S2, online only). As a sensitivity analysis (Table S3, online only), the cotrending test was also applied to groups consisting of (1) WMGHGs, TRF, and each , as well as to (2) TRF and each . In all cases, the conclusion is the same: at least part of the warming trend present in both regional and city‐level temperatures can be attributed to globally aggregated anthropogenic forcing.

Table 1 shows the estimated slope parameters B in Eq. (7) for each city. Most of the filtered temperature series () still contain a trend that in most cases indicates a faster warming than the regional rate (see also Fig. 2). These trends represent the excess warming relative to the large‐scale temperatures. This finding suggests that local factors had a considerable influence over the climate at the local scale and amplified the warming experienced regionally. The trends are significant for about 77% of the cities at the 10% level. The average slope coefficient over all cities shows an excess rate of warming of 0.44 °C (±0.12 °C SE) per century. A few cities experienced relative cooling up to −0.48 °C per century (Johnson City, Elkins, and Springfield), while most experienced higher warming than the regional average (e.g., New York, Allentown, and Milwaukee), with some having excess warming rates higher than 1 °C per century. While all weather stations selected by NOAA are intended to represent the climate conditions of cities, for some stations, the surroundings could be considered rural. Determining whether a weather station can be considered rural or urban is not trivial, and here we do not attempt such a formal classification. However, using Figures S3–S32 (online only), we identify eight stations with surrounding areas that are predominantly rural (Table 1). The average trend at these “rural” stations is 0.01 °C per century, while for the remaining stations, this value is 0.60 °C. This result can be interpreted as additional evidence supporting the fact that the proposed methodology is able to adequately separate global and local trends and to capture at least part of the urbanization effect. The presence of water bodies in the vicinity of weather stations can also influence the magnitude of the UHI. , In our sample, for most cities (23 out of 30), this measure is likely influenced by the presence of the sea, lakes, or rivers. However, these cities also show high levels of urbanization and for the vast majority, weather stations are located in highly urbanized surroundings. By contrast, only seven cities have no water bodies in their vicinities and, among those, only two are in urbanized areas. The average warming trend of stations in highly urbanized cities with water bodies is 0.56 °C, while this value is 0.07 °C in the case of less urbanized cities with no water bodies. Such results suggest that, in our sample, the urbanization effect dominates that induced by the presence of water bodies and that the combination of these two factors does not allow distinguishing the decrease in UHI commonly associated with the presence of water bodies.

Figure 2

Filtered local temperatures for the cities analyzed. Filtering of large‐scale trend and cycle was done by subtracting regional from local temperatures, as described by Eq. 6.

As discussed in the literature, the estimates of UHI can be more precise when stations are grouped according to the cities’ populations. , We consider groups of cities on the basis of population peak being larger/smaller than 200,000 inhabitants, yielding two subsamples of 15 cities each (Table 1). Large cities have an average excess warming rate of 0.56 °C (±0.12 °C SE) per century, while for smaller cities, this figure is 0.32 °C (±0.12 °C SE). The range for large cities basically excludes relative cooling (the lowest estimate being −0.06 °C per century), while for smaller cities, the lowest value is −0.48 °C in 100 years. The largest value for larger cities is 1.71 °C, while for smaller cities, it is 1.35 °C.

We emphasize that cannot be interpreted as a pure representation of the response to the UHI effect, but of a blend of multiple factors. The unobserved warming trend due to the UHI effect is distorted by , which may include the effect of factors other than UHI. However, as discussed below, correlation analyses between population and trend slopes and other characteristics suggest an association between UHI and the residual warming. The estimated magnitude of warming caused by local factors is similar to those reported in other studies that address the attribution of the UHI effect on temperature trends of urban stations. , , , The average rate of warming of the 50 most populous metropolitan regions in the United States was estimated to be 0.16 °C per decade over the period 1961–2010, with about 14 of them showing no warming or a slight decline. In the case of mainland China, the UHI effect at 45 stations was estimated in the range of 0.05–0.11 °C per decade, while for Beijing and Hebei, these estimates are 0.16 and 0.14, respectively.

Climate change is expected to modify the UHI, and some studies have projected that its intensity may decrease under very high emissions scenarios, but should remain similar to current days for moderate and low emissions scenarios (RCP4.5 and RCP2.6). A recent study based on temperature records for the period 2000–2015 proposed that the decrease in UHI warming caused by global climate change is already detectable in U.S. cities. For a longer sample (1985–2015), the authors report that they find no statistical differences between rural and urban stations. This result contrasts with what is found here and in much prior work, which has consistently found a positive UHI effect in observed air temperatures in the United States , , , and other parts of the world, , , , as well as with climate model projections in which significant decreases in the UHI intensity occur only under very high emissions scenarios for the end of the present century. The decrease in the UHI effect found by Scott et al. can be the result of low‐frequency oscillations that could dominate an underlying local warming trend in periods of time as short as 2000–2015. Moreover, for a given region, depending on the selection of rural and urban stations, different results can be obtained. , , ,

As mentioned in the introduction, urbanization is a dominant factor impacting the city's overall climate, and population count is commonly positively correlated with the UHI intensity. , , To provide insights about the origin of the excess warming rate observed in local temperatures, we calculated the correlation between the estimates of the coefficient B and the natural logarithm of population counts for each city, measured by the population peak (see Fig. 1 and Table 1, which also indicates when the peak occurred). As expected from the literature, the highest excess warming rates due to local factors occur in cities with the largest populations. The sample correlation (henceforth denoted by ρ) between local warming rate estimates and peak populations is 0.48 when considering all cities and 0.53 when only the statistically significant slope coefficients are included. Moreover, this relationship becomes stronger when only large cities are considered , while it is weaker for smaller cities . These results are depicted in Figure 3 and suggest that warming rates are associated with the extension and density of urban areas, as they are highly correlated with population counts. It is important to underline that these estimates suggest association, not causality. Moreover, the location of the stations may bias downward the estimated correlation coefficients, as some of them are in the outskirts of the cities where the UHI effect would be smaller. This could also explain why the correlation between the coefficient B and population is higher for large cities, as in most of them, the urban sprawl has integrated airports that were initially located in the outskirts. Correlation coefficients between average population growth rate and slope coefficients lead to a similar conclusion, but the magnitude of the correlations is lower, probably because the absolute value of population is better suited for inferring the UHI effect than using the growth rates.

Figure 3

Scatterplots of excess warming rates and peak population. Panel A shows the results for all cities, while panels B and C show the results for large cities (>200,000 inhabitants) and small cities (≤200,000 inhabitants), respectively.

Such relationships can be modulated by several other natural and anthropogenic factors represented by , which include city morphology, vegetation, elevation, where the station is located (i.e., close to downtown or in the outskirts), among others. In the case of smaller cities, the factors included in tend to lead to a higher variability in the estimated warming rates and thus to weaker correlations between population and urban warming (Fig. 3B and C). The differences in the magnitude of the correlation coefficients suggest that as cities become larger and more populated, the local effects tend to be dominated by the UHI. A natural determinant of the observed local warming rates is the elevation at which the station is located. The correlation between the estimated local warming trends and elevation is –0.54 when considering all cities, indicating that UHI effects are smaller as elevation increases. However, this relationship is much weaker for large cities than for small ones (with ρ = –0.64). These results reinforce the conclusion that as cities become larger, other factors that influence urban warming rates become less important.

The proposed methodology and results discussed above provide information about the effects of local factors on the urban warming trend through time. Features, such as level shifts created by local factors that modify the shape of the common warming trend, can provide additional insights to identify some of the most important drivers of the warming at the local scale. This can help to discriminate what factors in are more closely associated with the characteristics of than with those of . In some cases, distinguishing between a linear trend and a sudden or gradual shift in mean is difficult and both approaches can provide relevant information. The filtered local temperatures for each city depicted in Figure 2 suggest that some share similar features about the local warming evolution. The most notable feature is a possible level shift toward warmer temperatures during the mid‐20th century. To formally test for the existence of such level shifts in filtered temperatures, we apply the Bai and Perron methodology to test for multiple structural changes. Note that the representation of these breaks in the mean value as a sudden change is a convenient simplification. As can be seen from Figure 2, most of these changes are not abrupt but, in most cases, take more than a decade to attain their new level. Note that the tests will be consistent even in the presence of smooth or gradual changes. The estimates of the break dates are then interpreted as the date, within the period of change, that is most representative of the separation between the pre and postregimes. We applied this methodology to all city temperature series using regression (8). We used the AIC criterion to compare the linear trend and the shift in mean models and decide which one provides a better fit. The results in Table 1 show that nine of the analyzed temperature series contain a level shift toward warmer temperatures occurring around the mid‐20th century that is significant at the 5% level (using the Bai and Perron methodology). The magnitude of these shifts in mean is strongly correlated () with the estimated B coefficients. This reinforces the notion put forward that warming amplification and rapid changes in local temperatures are affected by common local drivers.

For 57% of the total cases, the shift in mean model produced lower AIC values than the trend model (Table 1). For most, a single shift was significant; only in three cases (Boston, Allentown, and Washington) was a second shift deemed significant. The average magnitude of the first mean shift is 0.33 °C. The estimated dates of the shifts in mean for the various cities are positively correlated with the dates when the population peak was reached (). In the case of cities having reached a peak in population (18 out of 30), the average peak date is 1957, while for cities for which a shift in the mean of was found (17 out of 30), the average break date is 1969. Many of the cities that show a shift in the mean of (11 out of 17) are smaller ones with peak populations lower than 200,000 inhabitants. For cities that show both a peak in population and a shift in mean, the correlation between the estimated dates for the shift in mean and the population peak is positive and very high (), although the number of observations is only 5. These results suggest that urbanization and urban sprawling would significantly slow down (or stop) some years after population count stabilizes and local warming would achieve a new level after the peak.

This may be clearer in smaller cities because the population peaks occur over a shorter span and are thus easier to detect given the time period available. Our results suggest that while population (and urbanization) is growing, the main effect of this growth is better captured by the estimated B coefficients. However, when population (urbanization) stabilizes, the effects of the level of urbanization attained are better expressed as a change in the mean. The overall results in this section strongly suggest that there is a strong association between the local warming, population count, and features that are city‐specific, such as whether a peak in population occurred.

Conclusions

We presented a methodology to conduct attribution studies for an urban temperature series that allows separating the contribution of global and local forcing factors. Cotrending tests are used to establish the existence of a common trend between globally aggregated TRF, its main global anthropogenic component, and large‐scale and local temperature series. We show that at least part of the observed warming trend in all cities has an anthropogenic origin. In the proposed methodology, the trend that corresponds to globally aggregated forcing factors is removed from the local temperature series and the existence of a local excess warming trend is evaluated and its features are characterized. Our results indicate that the vast majority of the cities (26 out of 30) analyzed show additional warming due to local factors and suggest that one of the main drivers is urbanization for which population count acts as a useful proxy. Increases in temperature due to urbanization and other local factors tend to level off after a city reaches its population peak, while in cities where population growth was uninterupted, temperatures maintain a trending behavior beyond the one attributed to global anthropogenic factors. Larger cities show a stronger relationship between local warming rates and population than smaller cities in which other factors, such as city elevation, seem to add more variability to local warming estimates.

Cities are important contributors to global greenhouse emissions and can also be particularly vulnerable to climate change. Given these characteristics and their exposure with reference to population and wealth concentrations, city‐level climate policy will have impacts not only at the local and regional levels, but also at the national and global scales. Understanding the different local and global drivers behind urban warming and the heterogeneity in the trends they generate is important to help decision makers design better mitigation and adaptation portfolios for risk reduction and management. Given the wide heterogeneity in the relation between regional and municipal warming, one could identify so‐called “green cities” as having a negative slope in the filtered trend and more so if the change in mean is negative; that is, a decreasing rate of excess warming even controlling for a positive increase in population. However, this is complicated by the common practice of placing weather stations in airports near the outskirts of the main urbanization center. This means that, while most cities in our analysis show additional local warming, the results obtained could be underestimating the UHI effect that corresponds to the central part of the cities. The problem is illustrated most notably by the case of Springfield (Illinois; Fig. S19, online only) and others, such as Lexington (Kentucky; Fig. S11, online only) in which the “green city” effect could be mainly attributed to the distance of the station to the city center. While the sample size is too small and a larger number of stations within a given city would be needed, it would be interesting to analyze in future work the features that may cause such differences with the majority of cities that experienced higher warming rate relative to the regional level. Also, the UHI effect shows strong seasonal and diurnal variation, with nighttime warming being commonly greater. , This is of particular importance since most of the local adaptation options to mitigate the UHI and the combined negative effects of global and local warming are less effective for the nighttime UHI, which is related to increasing health risks. Thus, it would be useful to apply the methods presented here to different seasons and to daytime and nighttime temperatures to provide a more complete characterization of the UHI effect.

Author contributions

F.E. and P.P. designed the study, analyzed the data, and wrote the paper. Both authors contributed equally to the study.

Competing interests

The authors declare no competing interests.

Figure S1. Regional area covered and location of weather stations in this study.

Figure S2. Temperature and radiative forcing time series.

Figures S3–S32. Location of the station in each city examined in this study.

Figure S33. Residuals, actual and fitted series from regressing on AMO.

Table S1. Location and details of weather stations used in this study.

Table S2. Bierens test for nonlinear cotrending around a linear trend for radiative forcing variables, regional and city‐level temperature series.

Table S3. Bierens test for nonlinear cotrending around a linear trend for radiative forcing variables and city‐level temperature series.

Click here for additional data file.