Urban heat islands in China enhanced by haze pollution

The urban heat island (UHI), the phenomenon of higher temperatures in urban land than the surrounding rural land, is commonly attributed to changes in biophysical properties of the land surface associated with urbanization. Here we provide evidence for a long-held hypothesis that the biogeochemical effect of urban aerosol or haze pollution is also a contributor to the UHI. Our results are based on satellite observations and urban climate model calculations. We find that a significant factor controlling the nighttime surface UHI across China is the urban–rural difference in the haze pollution level. The average haze contribution to the nighttime surface UHI is 0.7±0.3 K (mean±1 s.e.) for semi-arid cities, which is stronger than that in the humid climate due to a stronger longwave radiative forcing of coarser aerosols. Mitigation of haze pollution therefore provides a co-benefit of reducing heat stress on urban residents.

T he urban heat island (UHI) represents one of the most pronounced surface climate changes caused by human activities 1 . The mechanism underlying the UHI formation is generally thought to be biophysical in nature, arising from large differences between rural and urban land in surface properties, including sensible heat dissipation or convection efficiency, evaporative cooling, sunlight reflection and artificial heating 2,3 . The main contributors to the daytime UHI are reductions in sensible heat convection efficiency and evaporative cooling in urban land [3][4][5] . At night, the heat released from energy use and from solar energy stored in buildings is a major biophysical factor responsible for urban warming. The UHI increases the number and the intensity of heat waves in cities, thus aggravating heat stress on urban residents 6 .
Cities are also the dominant source of anthropogenic aerosols having large impacts on the biogeochemistry of the atmosphere 7,8 . The haze plume formed from urban aerosols alters regional precipitation patterns outside the city 9 and contributes to radiative forcing on the global climate 10 . The haze biogeochemical effect on the climate of urban land itself is, however, still not well understood, in large part because of the difficulty in disentangling the opposing effects aerosols have on the surface shortwave and longwave radiation 11 . Aerosols generally reduce the amount of shortwave radiation reaching the ground surface, creating a cooling effect on the surface. However, they are much more effective in absorbing and emitting radiation than water vapour and greenhouse gases in the longwave atmospheric window (wavelength range 8-12 mm) under specific conditions, thus having the potential to increase the longwave radiation energy received at the urban land surface 12 . The overall effect depends on the initial particle size and size growth due to ageing and absorption of water vapour 13 .
So far, a quantitative evaluation of the haze contribution to the UHI has not been attempted. There are two difficulties. First, standard urban land surface models 14 do not include parameterizations for incoming shortwave and longwave radiation, because these are properties of the atmosphere aloft, not biophysical properties of the surface itself. Second, the aerosol radiative forcing is not a prognostic variable in atmospheric data assimilation models and in most climate models.
Cities in China are burdened with unprecedented levels of aerosol pollution. Here we present an empirical analysis using satellite observations to show that urban haze pollution is a contributor that intensifies the UHI in China at night. In the analysis, the surface UHI intensity DT is the difference in surface temperature between the urban and the adjacent rural land 3 , and haze pollution is measured by the aerosol optical depth (AOD). We then use an urban climate model in conjunction with an observation minus reanalysis (OMR) method for aerosol longwave radiation enhancement, to quantify the haze contribution to the surface UHI intensity in three climate zones (humid, semi-humid and semi-arid) across China. We find that haze pollution intensifies the nighttime UHI in China through an increase in incoming longwave radiation. This warming effect is greater in semi-arid cities compared with cities in humid and semi-humid climates.

Results
Drivers of UHI spatial variations. We analysed the annual mean DT measured by the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument on board of the Aqua satellite from 2003 to 2013 for 39 cities across Mainland China (Fig. 1a). Being a proxy for heat release from anthropogenic sources and from solar energy stored in buildings, population is a predictor frequently used to explain city-to-city variations in the nighttime UHI intensity observed by satellites 4,15-17 and by weather stations 2 . Our results indicate that, in contrast to studies reported for other regions of the world 4,15 , city population is actually a poor predictor of the nighttime DT variations among these cities (Fig. 1b). Shanghai, located in southeast China and the largest city we analysed (population 14 million), shows a weak surface UHI (1.5 K), whereas Hami, a small city (population 0.45 million) in northern China, exhibits one of the strongest surface UHIs (5.0 K). The overall correlation between population and nighttime DT is not statistically significant (P ¼ 0.17).
Another unusual feature is the diurnal variations of the surface UHI. The MODIS-derived nighttime DT (3.4±0.2 K, mean±1 s.e.) is higher than the daytime value (2.1 ± 0.3 K; Po0.001). The diurnal contrast is especially striking for cities in the semi-arid climate, where the mean nighttime DT is 4.0 ± 0.4 K, but the daytime DT is only 0.3±0.5 K (Po0.001; Fig. 2c,f). The diurnal patterns in China differ from those observed by satellites for North America 18 , Europe 19,20 , South America 20 and Oceania 20 where the daytime surface UHI is stronger than the nighttime UHI and where the UHI of semi-arid cities is generally weak at night 4,17,18 . Our results are broadly consistent with the UHI spatial pattern documented for China in a previous study using a shorter MODIS time series 21 , although our UHI intensity is generally greater, because we only used pure urban and rural pixels to calculate DT.
One explanation for these unique surface UHI patterns is related to haze pollution. We find that the spatial variations of the annual mean nighttime DT are significantly correlated with the difference in AOD between urban areas and the adjacent rural land ( Fig. 1c; DAOD, urban AOD minus rural AOD; Po0.01). Significantly positive correlation is also found between the summer nighttime DT and DAOD (Po0.01). Cities having a thicker haze layer than the surrounding rural environment tend to display a stronger UHI. We use DAOD, because DT is also a perturbation signal in reference to the rural background. The AOD itself is not a good predictor of the DT variations (P ¼ 0.48). Only after controlling for DAOD does DT show significant dependence on population (Po0.01).
There is no evidence of haze enhancement on the daytime DT. The correlation between the annual daytime DT and DAOD is poor (P ¼ 0.43; Fig. 3). Repeating the correlation analysis for the summer season reveals similarly poor correlation (P ¼ 0.50). Instead, the most important factors explaining the daytime DT variations are population, urban-rural difference in normalized difference in vegetation index and cropland fraction of the rural background (P-valueso0.001). Annual mean precipitation exerts a strong control on the daytime DT in North America (Po0.001; ref. 3) but a weak control in China (P ¼ 0.06). This regional difference can be explained by the fact that cropland is a more prominent non-urban land cover in China than in North America. Irrigation is commonplace in China, with 48% of the farmland receiving water from irrigation in addition to water supplied by rain 22 . Domesticated plants supported by irrigation water do not behave in the same way as natural ecosystems in terms of surface energy exchanges. After excluding cities whose adjacent rural area consists of 450% cropland pixels, precipitation becomes a significant controlling factor (linear correlation ¼ 0.57, Po0.01, number of observations ¼ 21). The UHI dependence on precipitation and irrigation highlights the important role of the rural background environment in calculating DT.
Attribution of the haze effect. Figure 1c can be viewed as empirical evidence supporting the long-held hypothesis that urban haze pollution is a contributor that intensifies the UHI 2 . We now make a quantitative attribution of the haze contribution to the nighttime DT by combining climate model calculations with analysis of surface longwave radiation observations. The surface radiation data are used in conjunction with the surface longwave radiation calculated by an atmospheric data assimilation model, to obtain the sensitivity of L k to AOD, that is, the amount of enhancement in L k in response to a unit increase in AOD (Methods). Here, L k is the downward longwave radiation received by the surface including emissions and scattering of air molecules and aerosols. In the climate model, the urban land is parameterized as a separate land unit at the subgrid level. We force the model with an assimilated atmosphere and save the surface energy balance variables of urban and non-urban subgrid land units for offline UHI attribution 4 . The attribution method separates the contributions of external radiative forcing, energy redistribution via aerodynamic resistance-associated sensible heat convection and energy redistribution via evaporation 23 .
In this framework, the aerosol effect is an external forcing similar to anthropogenic heat release and to changes in the surface shortwave radiation arising from the urban-rural surface albedo difference, and can be expressed as, where (DT) h is haze contribution to the UHI intensity, l 0 (E0.20 K m 2 W À 1 ) is the local intrinsic climate sensitivity, f is a dimensionless energy redistribution factor and DL k is the urban-rural contrast in L k calculated as the product of the satellite-observed DAOD and the longwave radiation sensitivity to AOD. Our AOD sensitivity values (Table 1) fall in the range of those calculated with radiative transfer models [24][25][26][27] . According to the observations of a ground-based aerosol remotesensing network 28 , the aerosol Ångström exponent is smaller in the semi-arid northwest Chinese cities, indicating larger particle sizes, than in cities in the humid central and eastern China. The sensitivity for the semi-arid climate zone is much higher than for the humid climate zone, confirming a stronger longwave radiative forcing of coarser particles 11,25 . Our estimates of DL k , B1.1 and 8.0 W m À 2 for the cities in the humid and semi-arid climate zone, respectively, are lower than those reported from paired observations at urban and rural sites 29,30 , because we did not consider the L k enhancement caused by a warmer urban boundary layer 31 and emissions from urban canopy walls 32 . In the model domain, L k represents the downward longwave radiation incident on a reference plane above the urban canopy, which is the first model grid height. We estimate that the haze contribution to the nighttime DT is 0.70 ± 0.26 K (mean ± 1 s.e.) for the semi-arid cities and is small for the other two groups of cities ( Table 1). The larger (DT) h in semi-arid climate is a result of less efficient energy redistribution (smaller f-values due to larger aerodynamic resistance and Bowen ratio; equation (5)) between the land and the atmosphere, a larger urban-rural contrast in pollution level and a stronger longwave radiative forcing of coarser aerosols.
In the semi-arid climate zone in China, both the urban and its adjacent rural area are affected by coarse mineral particles transported from the Taklimakan Desert and the Gobi Desert 33 . In the urban environment, road fugitive dust, constructionderived dust, and domestic heating and cooking are additional sources of coarse mode aerosols, explaining why the ratio of PM 10 to PM 2.5 concentration is greater in the semi-arid cities in Northwest China than in the humid cities in South China [34][35][36] . The urban-rural AOD difference observed here appears to result from these urban anthropogenic sources.
The haze contribution to the daytime DT is uncertain because of the opposing effects aerosols have on the surface shortwave and longwave radiation 11 , but the lack of correlation between the daytime DT and DAOD (Fig. 3) suggests that it may be negligible. This inference is consistent with the model results. The daytime DT determined online by the Community Land Model (CLM) model, denoted as CLM in Fig. 2, is in good agreement with the MODIS-derived values (denoted as MODIS) for the cities in humid and semi-humid climates (Fig. 2a,b). For the cities in semi-arid climate, the model online result overestimates the observation, but the offline diagnostic calculation (denoted as Calculated), which is the sum of all the terms in equation (2), shows a good agreement (Fig. 2c). The overall agreement leaves little room for an additional contribution due to the haze effect, implying that the relative reduction of the shortwave radiation in the city in reference to the rural background is roughly equal to the relative enhancement of the longwave radiation. This offsetting effect of aerosols on radiation has been reported previously by ref. 37 and the year-long observations at an urbanrural site pair support this interpretation 29 . Comparison of the model results and the MODIS observation for the summer also reveals very good agreement for the daytime.
According to the attribution diagnostics, the main contributor to the daytime UHI in the humid climate is the reduction in sensible heat convection efficiency of the urban land, not a reduction in evaporation. In the semi-arid climate, the role of convection is reversed, contributing to a cooling signal. These results are in agreement with those obtained previously for cities in North America 4 .
In contrast to the daytime results, the modelled nighttime DT is too low in comparison with the MODIS observations ( Fig. 2d-f), even though the same model is able to reproduce the observations in North America. The atmosphere in the model domain is free of haze pollution. With the inclusion of the haze contributions calculated offline (Table 1), the modeled DT is still biased low. One reason for the low bias is that the model scheme does not have a complete accounting of all sources of anthropogenic heat release in the urban land 38 . The anthropogenic heat flux is an important contributor to the nighttime UHI 21 . Another possibility is that equation (1) has omitted a dynamic mechanism  associated with the haze pollution. Aerosols are known to warm the atmosphere 12,39 , potentially making the boundary layer air above the urban land more stable than that above the rural land and thus reducing the efficiency of energy dissipation from the urban surface to the atmosphere. The end result is an amplified UHI intensity. It appears that this stability mechanism is especially strong in the humid climate, as suggested by the large model bias error (Fig. 2d), although a definite answer will require an improved anthropogenic heat parameterization for China.
Our study implies that abatement of haze pollution has a co-benefit of reducing heat stress on urban residents. The UHI intensity is a perturbation signal in reference to the rural background temperature. A complete assessment of the haze effect must recognize that the rural atmosphere is also changing. The MODIS observation indicates that the rural AOD in China is 0.20 and 0.53, greater than that in North America in the semi-arid and the humid climate zones, respectively. The data in Table 1 suggest that the rural land in China may be receiving B15 W m À 2 more longwave radiation energy at night than under haze-free conditions. Our subgrid energy balance analysis is not suited for quantifying the impact of the rural background change, because the change signal is regional in scale. The nighttime temperature in China increased at 0.47 K per decade from 1979 to 2012 (ref. 40), which is roughly twice the global mean temperature trend for the same time period 41 . We hypothesize that haze pollution has contributed to the accelerated warming in China. If the hypothesis is proved valid, pollution abatement should also relieve heat stress on rural populations.

Methods
Satellite data. We selected 39 cities in mainland China with the consideration that there is at least one city for each province (Supplementary Table 1). More cities were added for the three largest provinces (Xinjiang, Qinghai and Inner Mongolia), to ensure even distribution across the country. In our selection, cities located in mountainous regions were avoided. Of the selected cities, 19 are in humid climate, 9 in semi-humid climate and 11 in semi-arid climate, according to the Köppen-Geiger climate classification 42 .
We used the (MODIS) Aqua land surface temperature (LST) product (MYD11A2), to retrieve paired urban and rural LST. The MYD11A2 product comprises 8-day clear-sky LST observations at 13:30 and 1:30 local time, with a 1 km spatial resolution. The study period is from 2003 to 2013. A set of 3 Â 3 pixels was chosen from the urban core, except for three small cities for which we were only able to obtain one to three pure urban pixels. The surrounding rural areas were represented by up to four 3 Â 3 pixel patches at the four sides of the city, but mountain and water pixels were avoided. The selected urban-rural pixel pairs do not differ by, on average, 4100 m in elevation and by 40.1°in latitude.
Three features of the MYD11A2 product make it particularly suitable for the UHI detection. First, a MODIS cloud mask (MYD35) has been applied to filter out cloudy conditions; thus, cloud interferences are avoided. Second, a generalized split-window algorithm using two longwave bands in the atmospheric window is used to correct atmospheric water vapour and haze effects, and to reduce the sensitivity to errors in the surface emissivity 43 . The average bias (MODIS LST minus ground-based measurement of LST) is À 0.15 ± 0.73 K in conditions of low haze pollution 44 (mean±1 s.d., number of sites n ¼ 6) and is not different statistically (two-sample t-test, P ¼ 0.98) from the mean bias of À 0.17 ± 1.66 K (n ¼ 4) in conditions of high haze pollution 45 . Third, the brightness temperature has been corrected for surface emissivity to obtain the true LST. The surface emissivity data, developed by a MODIS science team 46 , delineate the land surface into 17 categories (including a category for urban) and account for the bidirectional radiation distribution factor of each category. The impact of emissivity on the surface temperature retrieval has been discussed elsewhere 43,46 .
Land cover change over time can complicate satellite UHI observations 47 . Steps were taken to ensure that the selected urban and rural pixels stayed as urban and rural throughout the duration of the study. The selected pixels were first validated by the MODIS land cover product (MCD12Q1, resolution 500 m) with the IGBP classification for year 2010 and then cross-checked against Google Earth. As no abandonment of urban land has occurred in these Chinese cities, these rural pixels should be rural in the years before 2010. The urban pixels were from the urban core; hence, they should be urban in the prior years as well as in the future years. To address the question of whether the pixels designated as rural in 2010 were converted to urban in the later years, we compared our selected rural pixels for each city against the MODIS classification for year 2013, the last year of our study period. Only the rural pixels for eight cities experienced some change from 2010 to 2013, but the change is small: o5% of the rural pixel selected based on the 2010 classification was converted to the urban class in 2013.
The MODIS Level 2 aerosol product (MYD04_3K) was used in this study. Compared with the nominal 10 km aerosol product, this product has a finer resolution of 3 km and is more suitable for determining urban and rural aerosol characteristics contrast. The AOD values were determined with the dark target algorithm 48 . Only high-quality (quality flag ¼ 3) data were used. For each pair of city and its adjacent rural area, we calculated the averaged AOD from 2003 to 2013 and subtracted the rural mean AOD from that of the urban area to obtain DAOD. The community land model and model simulations. The model we used is the CLM (version 4.0) of the NCAR's Community Earth System Model (CESM) system 49 . Each gridcell in the CLM is configured with five land units: glacier, wetland, vegetated, lake and urban 50 . In this study, the urban and vegetated land units in the same gridcell were used to represent an urban-rural site pair. The model was driven by the atmospheric forcing data described in ref. 51. The simulation period  Offline attribution of UHIs. We adopted the attribution method first developed for the investigation of local deforestation effects 23 and later extended to the study of the surface UHI of North American cities 4 . Briefly, the CLM subgrid model outputs were analysed in a surface energy balance framework that separates the contributions to the urban-rural temperature difference according to five biophysical drivers: net shortwave radiation or surface albedo change, change in sensible heat dissipation efficiency, evaporation reduction, heat storage change and anthropogenic heat release. Their contributions to DT are expressed as: where f is an energy redistribution factor, which is a function of Bowen ratio and the aerodynamic resistance (equation (3) below), l 0 is a local climate sensitivity parameter (l 0 ¼ 1/4sT 3 E0.2 K W À 1 m 2 , where s is the Stephan-Boltzmann constant) and the terms on the right side of the equation represent contribution associated with difference in the net apparent radiation (term 1), sensible heat convection efficiency (term 2), evaporation (term 3), heat storage (term 4) and anthropogenic heat flux (term 5) between the urban and the rural area of the same grid. The change terms Df 1 and Df 2 are expressed as: The surface UHI intensity was computed in two different ways. First, DT was calculated online by CLM as the surface temperature difference between the urban and the vegetated land unit in the same gridcell (the green bars in Fig. 2). In the second method, DT was computed as the sum of the individual contributions in the offline analysis according to equation (2) (the blue bars in Fig. 2). The atmosphere in CESM is free of haze pollution, but the pollution effect can be quantified using the offline diagnostics according to equation (1). The most relevant diagnostic variable is the energy redistribution factor, f, which describes the efficiency of energy redistribution between the surface and the atmosphere, and is defined as, where r is air density, C p is specific heat of air at constant pressure, r a is aerodynamic resistance to sensible heat, s is the Stephan-Boltzmann constant, T s is the surface temperature and b is Bowen ratio calculated as the ratio of the subgrid sensible heat to latent heat flux. These diagnostic variables are provided by the CLM model for every model grid. In a model grid where the f-value is higher, the surface UHI intensity will probably be lower, because energy is dissipated more efficiently from the surface to the atmosphere.
To estimate the haze enhancement on the surface UHI, we need to quantify the difference in the incoming longwave radiation between the urban and the rural area (DL k ). We first determined the sensitivity of L k to AOD, that is, the amount of enhancement in L k in response to a unit increase in AOD, using the OMR method 52 . The reanalysis longwave radiation data were provided by the Modern-Era Retrospective Analysis for Research and Application. The observational data came from three ground sites in the ChinaFLUX network 53 representing the three climate zones (Yongfeng, 32.21°N, 118.67°E; Luancheng, 37.88°N, 114.68°E; Xilinguole 43.53°N, 116.66°E). Yongfeng (Jiangsu Province) and Luancheng (Hebei Province) are cropland sites and Xilinguole is a grassland site in Inner Mongolia. As Modern-Era Retrospective Analysis for Research and Application does not consider haze pollution but the observation was impacted by haze, some of the difference between the observed and reanalysed L k can be considered to result from haze pollution. In this study, the regression of the OMR L k against the MODIS AOD is given as [OMR L k ] ¼ a þ b Â [AOD], where the intercept a represents the reanalysis model bias 54 and the slope parameter b was taken as the sensitivity of L k to AOD (Table 1). Next, the urban L k enhancement was calculated as the product of the AOD sensitivity and the observed AOD difference between the paired urban-rural pixels. Finally, the haze contribution to the UHI was calculated according to equation (1). The s.e. of f and DL k were calculated after averaging the values for cities belonging to the same climate zone. Uncertainty on (DT) h was determined with a Monte Carlo procedure using 100,000 realizations and a Gaussian error distribution for all input variables.
Data availability. All the satellite data used in this study can be downloaded on the MODIS product website (https://lpdaac.usgs.gov/dataset_discovery/modis/ modis_products_table/myd11a2). Other relevant data in this study are available from the authors.