-
The Community Land Model version 5.0 (CLM5.0) is the land component of the Community Earth System Model (CESM) (Danabasoglu et al., 2020). It uses a nested hierarchy consisting of up to five land units (glacier, urban, agricultural, vegetation, and lake) in a grid cell to represent the land surface heterogeneity at the subgrid level (Lawrence et al., 2019a). The vegetation unit is further broken down to plant functional types. In the present study, these plant functional types are regrouped to tree, grass, shrub, and bare soil tiles. The surface fluxes are calculated for each tile at subgrid level, then the area-weighted average is computed at grid level. A full description of the model data is provided by Zhang et al. (2023). A brief summary is given below.
The simulation was conducted at a 0.9°×1.25° resolution under the Representative Concentration Pathway 8.5 (RCP8.5) scenario from 2015 to 2100. There are 355 grid cells containing both urban and rural tiles across China. They fall into three climate zones: humid (201 grids), semi-humid (118 grids), and semi-arid (36 grids) (Fig. 1) according to the Köppen–Geiger climate classification. In CLM5.0, plant phenology is prescribed by satellite observations (Zhang et al., 2023). The surface data were fixed at present-day levels in the simulation, which used the urban land cover in 2020 under the Shared Socioeconomic Pathway 5 (SSP5) scenario from He et al. (2021). Prior to the coupled simulation, a spin-up simulation was conducted offline using GSWP3 (Global Soil Wetness Project phase 3) observation data from 1991 to 2010. The forcing data were cycled for 60 years until the surface climates reached equilibrium [Fig. S3 of Zhang et al. (2023)]. The initial land condition for the spin-up simulation is a default initial file in the year 2011 provided by CESM2. Model outputs are archived for eight land tiles (urban, rural, tree, grass, shrub, bare soil, crop, and lake), including key variables on the physical state, surface energy fluxes, and atmospheric forcing conditions. The urban tile is based on the urban canyon concept, which consists of five components: roof, sunlit wall, shaded wall, pervious canyon floor, and impervious canyon floor. Hourly outputs are available for two periods: 2019–23 and 2096–2100. We only used the variables of urban and rural subgrid tiles at hourly intervals from 2019 to 2023. Unless stated otherwise, our analysis is restricted to midday [1300 LST (local solar time)] and midnight (0100 LST). The two selected times can represent the typical daytime and nighttime conditions according to the AUHI diurnal variations [Fig. S1 in the electronic supplementary material (ESM)]. The broad spatial patterns among the three climate zones are unchanged if longer averaging periods (1200–1600 LST and 0000–0400 LST) are used (Table 1). Also archived are subgrid auxiliary data, including LAI, urban street height-to-width ratio, urban area fraction, and plant functional types.
Averaging period (LST) Humid Semi-humid Semi-arid Daytime 1200–1300 0.34 ± 0.14 0.50 ± 0.14 0.73 ± 0.13 1200–1600 0.46 ± 0.12 0.58 ± 0.10 0.73 ± 0.13 Nighttime 0000–0100 0.75 ± 0.31 0.95 ± 0.31 0.81 ± 0.37 0000–0400 0.71 ± 0.30 0.91 ± 0.30 0.80 ± 0.36 Daily 24-h mean 0.58 ± 0.19 0.75 ± 0.17 0.72 ± 0.23 Table 1. Daytime, nighttime and daily mean AUHI (mean ± one standard deviation; K) for three climate zones across China.
The AUHI is computed as the difference in the screen-height air temperature between the urban tile and the rural tile in the same model grid. The temperature of the rural tile is an area-weighted average of bare soil, crop, tree, grass, and shrub temperatures. In Fig. S2 in the ESM, we show the air temperature of these individual tiles for three grid cells in three different climate zones, corresponding to Nanjing, Beijing and Urumqi, to demonstrate typical within-grid air temperature variations.
-
In our study, the AUHI intensity is quantified as the 2-m air temperature difference between the urban and rural tiles (∆T2) in the same model grid cell. In the following, we describe a diagnostic framework for attributing ∆T2 to different biophysical drivers. The basis of this framework is the analytical solution of surface temperature Ts from the surface energy balance equation (Zhao et al., 2014),
where Tb is air temperature at the blending height,
$R_n^*$ is apparent net radiation, QA is anthropogenic heat flux, Qs is storage heat flux, f is the dimensionless energy redistribution factor, and${\lambda _0}$ is the local climate sensitivity. Qs is calculated as the residual of the surface energy balance equation (Lawrence et al., 2019b). It is the sum of all the heat storage components, including biomass heat storage (Swenson et al., 2019), soil heat storage, and heat storage in canopy air. The anthropogenic heat in CLM5.0 consists of space heating and air conditioning loads to keep the indoor temperature within a comfortable range. Because this parameterization omits vehicle heat release, the anthropogenic heat flux in CLM5.0 is biased low (Zhang et al., 2023). Details of the parameterization are described by Oleson and Feddema (2020). The three intermediate variables are given bywhere σ is the Stefan–Boltzmann constant,
$\varepsilon $ is surface emissivity, K↓ is incoming solar radiation, L↓ is incoming longwave radiation, α is surface albedo, rt is the total resistance for heat transfer between the surface and the blending height, β is the Bowen ratio (ratio of sensible heat flux to latent heat flux), ρ is air density, and Cp is the specific heat of air at constant pressure.On the assumption that sensible heat flux is constant with height between the surface and the blending height, we can relate T2 to Ts as (Wang and Li, 2021)
where R = ra / rt, and ra is the aerodynamic resistance between the 2-m height and the blending height.
Differentiation of Eq. (5) yields a diagnostic equation for the AUHI:
where Δ denotes the urban–rural difference in a variable. The changes in surface roughness and Bowen ratio can cause changes in energy redistribution. Their individual contributions are given by
Terms on the right hand side of Eq. (6) represent contributions from the urban–rural difference in convection efficiency (term 1), evaporation (term 2), radiative forcing (through changes in surface albedo and emissivity, term 3), heat storage (term 4), and anthropogenic heat release (term 5).
In our diagnostic analysis, the blending height is equivalent to the first model grid height (~30 m). The two resistance terms are calculated from the following diagnostic relations:
where H is sensible heat flux. All other variables in Eq. (6) are provided by the model.
To validate the attribution framework, we compared the modeled AUHI with the sum of the component contributions calculated offline according to Eq. (6). The results (Fig. S3) show excellent agreement at midday (r = 0.99, p < 0.001, RMSE = 0.03 K) and acceptable agreement at midnight (r = 0.75, p < 0.001, RMSE = 0.26 K). Expressed as climate zone means, excellent agreement is found for midday, with the absolute errors less than 0.03 K. At midnight, the absolute errors are less than 0.21 K.
-
In this study, three statistical analysis methods were used to determine the main drivers of spatial patterns of the AUHI across China. First, we performed a single-variable correlation analysis to find the correlation coefficients between the midday and midnight AUHI and two types of drivers: (1) six background climate drivers [annual mean air temperature T, annual precipitation P, incoming solar radiation K↓, incoming longwave radiation L↓, net radiation Rn, and vapor pressure deficit (VPD)]; (2) a background ecological driver (rural LAI). All the drivers are provided by the CESM2 outputs. Danabasoglu et al. (2020) evaluated the performance of CEMS2 in simulating historical air temperature and precipitation.
Second, dominance analysis (Budescu, 1993; Azen and Budescu, 2003) was used to determine the relative contribution of each driver to the spatial variations in the AUHI. Prior to this analysis, all the variables were normalized between 0 and 1, with 0 corresponding to the minimum value and 1 to the maximum value. Firstly, a stepwise multi-variable regression was used to screen the predictors. Let us suppose that we have m predictors after this screening. The regression model with all m predictors is the complete model. We can establish other possible regression models with fewer than m predictors (including null). These models are called subset models. There are a total number of 2m − 1 possible subset models. To determine the dominance of predictor X1, we first select models without X1 as a predictor. We then calculate the increase in variance by adding X1 to the selected subset models. This process is repeated for other m − 1 predictors. The relative contribution of X1 is expressed as the ratio of variance increase by X1 to the total variance increase by all m predictors.
Third, a covariance analysis (Zhao et al., 2014) was performed to determine how the main driver obtained with the above two methods interacts with the biophysical contributions to influence the spatial variations in the AUHI. Let Cc, Ce, Cr, Ch and CAH represent contributions from the five biophysical contributions [terms 1 to 5 in Eq. (6)]. Equation (6) can be rewritten as
where e denotes a residual error resulting from nonlinear interactions among the biophysical factors. The spatial covariance of the AUHI and the main driver is equal to the sum of the covariance between each component contribution and the driver:
Averaging period (LST) | Humid | Semi-humid | Semi-arid | |
Daytime | 1200–1300 | 0.34 ± 0.14 | 0.50 ± 0.14 | 0.73 ± 0.13 |
1200–1600 | 0.46 ± 0.12 | 0.58 ± 0.10 | 0.73 ± 0.13 | |
Nighttime | 0000–0100 | 0.75 ± 0.31 | 0.95 ± 0.31 | 0.81 ± 0.37 |
0000–0400 | 0.71 ± 0.30 | 0.91 ± 0.30 | 0.80 ± 0.36 | |
Daily | 24-h mean | 0.58 ± 0.19 | 0.75 ± 0.17 | 0.72 ± 0.23 |