Seasonal and Diurnal Variations in the Relationships between Urban Form and the Urban Heat Island Effect

Abstract

At the city scale, the diurnal and seasonal variations in the relationship between urban form and the urban heat island effect remains poorly understood. To address this deficiency, we conducted an empirical study based on data from 150 cities in the Jing-Jin-Ji region of China from 2000 to 2015. The results derived from multiple regression models show that the effects of urban geometric complexity, elongation, and vegetation on urban heat island effect differ among different seasons and between day and night. The impacts of urban geometric factors and population density in summer, particularly those during the daytime, are significantly larger than those in winter. The influence of urban area and night light intensity is greater in winter than in summer and is greater during the day than at night. The effect of NDVI is greater in summer during the daytime. Urban vegetation is the factor with the greatest relative contribution during the daytime, and urban size is the dominant factor at night. Urban geometry is the secondary dominant factor in summer, although its contribution in winter is small. The relative contribution of urban geometry shows an upward trend at a decadal time scale, while that of vegetation decreases correspondingly. The results provide a valuable reference for top-level sustainable urban planning.

1. Introduction

Urban heat islands (UHIs) magnify the intensity of heatwaves globally [1,2] and impose a disproportionately large energy burden on cities [3]. The scientific design of urban form has been considered a key factor to address these issues due to its profound impacts on the urban thermal environment and its strong “lock-in” effect [4]. As a result, the urgency of urban form optimization in addressing global sustainability has been emphasized both by The 2030 Agenda for Sustainable Development and The New Urban Agenda [5,6].

Progress has been made in examining the relationship between urban form and the UHI effect at the neighborhood scale and city scale. At the neighborhood scale, both the local climate zone (LCZ) and precinct ventilation zone (PVZ) are key components of urban form that influence land surface temperature. For instance, based on the use of LCZ, Yang et al. showed that urban architectural patterns were one of the important drivers of climate change, and high-density high-rise buildings can increase surface temperatures [7]. Jie He, Lan Ding, and Deo Prasad developed the PVZ concept and conducted field measurements to verify the relationships among urban form indicators such as compactness, building height, street structure, wind patterns, the UHI effect, and outdoor thermal comfort [8]. However, the relationships among these factors have been found to be dependent on urban size, climate background, and, in particular, the study scale [9,10,11]. Since there may be significant differences in the relationship between urban form indicators (UFIs) and UHI effects at the neighborhood and city scales [12], researchers have called for additional approaches or cross-sectional empirical studies to help clarify the ambiguous effects of urban-scale characteristics [13,14]. City-scale studies have been conducted worldwide, but these studies either failed to take seasonal and diurnal variations into consideration or lacked sufficient consideration of city-scale urban form factors [9,15,16].

At the city scale, urban forms include urban shape geometric indicators (GEO), urban size and density indicators (SIZ), and urban vegetation indicators (VEG). Specifically, it includes shape continuity (CONTIG), elongation (ELONG), fractal dimension (FRAC), urban area size (AREA), urban development index (UDI), population density (POP), normalized differential vegetation index (NDVI), and tree cover (TREE). These urban form indicators not only directly influence UHIs by affecting the physical properties of the underlying surface but also indirectly influence UHIs by affecting traffic patterns, greenhouse gas (GHS) and pollution emissions, air circulation, and other socioeconomic processes [17,18,19,20,21,22,23,24,25]. As Figure 1 shows, urban forms influence UHIs in different ways by affecting various energy flows. Due to the distinct physical and biological characteristics of urban forms in different seasons and between day and night, the magnitude of these positive and negative effects will change in a complex manner at the city scale. According to the urban energy balance theory, all the aforementioned effects of urban form ultimately affect the partitioning of the energy balance, which in turn modifies the urban microclimate. An adopted energy balance equation is provided by Oke [26] as Equation (1), where all terms are flux densities and Q_r is the net all-wave radiative flux; Q_T is the anthropogenic heat (AH); Q_E is the convective or turbulent sensible heat flux; Q_L is the latent heat flux; Q_S is the stored energy; Q_A is the net horizontal heat advection. Due to the great seasonal and diurnal variations in the energy balance of urban surfaces [27], as evidenced by variations in the ratios of Q_r to Q_T and of Q_E to Q_A (Figure 1), the effects of UFIs on UHIs among different seasons and times of day will vary. Lack of understanding of the heterogeneity in the relationship between UFIs and the UHI effect on a macroscale will hinder the formulation of more targeted planning policies for sustainability.

Q_r + Q_T = Q_E + Q_L + Q_S +Q_A

(1)

In addition to the deficiencies discussed above, further gaps exist. One problem is the research scale [30,31]. Some researchers have claimed that urban size and density have little impact on the UHI effect in studies at the global or national scale [15,32] but have a great impact in studies at the regional or smaller scale with a consistent geographical background [33,34]. A scale that removes the nonurban effects is believed to be appropriate in the study of the urban form effect on UHIs, but this approach has not been widely adopted [35]. Another issue is associated with data structure and statistical models. Limited by the quantity of nonexperimental urban-scale data, it is difficult for a commonly adopted single-period, cross-sectional data model to determine causality (Wooldridge, 2014). In addition, the selection of urban form variables in the modeling of UHI intensity is rarely demonstrated reliably, which may lead to important bias. In general, univariate models are vulnerable to missing variables, while the multivariate models are challenged with the multiple collinearity or loss of specific variable coefficient estimates caused by the variable-discarding strategy.

The purpose of this paper is to study the seasonal and diurnal variances in the relationship between urban form and the UHI effect at a macroscopic city scale. Several improvements were made to address the aforementioned gaps. First, a region with a uniform climatic background and densely distributed cities (the Jing-Jin-Ji region of China) was chosen as our research area. Second, we extracted pooled cross-sectional data for 150 cities from 2000 to 2015 in four periods every five years for modeling to consider the temporal changes in the data [36]. Then, ridge regression (RR) models were constructed through parameter optimization to obtain a more reliable estimation of the impact of urban forms on UHIs. Finally, multi-perspective statistics, including the correlation coefficient, regression slope coefficient, and relative contribution rate, were applied to obtain a more detailed understanding of the seasonal and diurnal variations in the relationships between urban form and UHIs.

2. Materials and Methods

2.1. Study Area

Located in northern China, the Jing-Jin-Ji region is one of the largest urban agglomerations and the most developed regions in China (Figure 2a). The Jing-Jin-Ji urban agglomeration includes the cities of Beijing, Tianjin, and Hebei Province, which includes Shijiazhuang, Tangshan, Qinhuangdao, Handan, Xingtai, Baoding, Zhangjiakou, Chengde, Cangzhou, Langfang, and Hengshui, all of which have different development characteristics. The region covers 217,156 km², 1.9% of China’s territory. The Jing-Jin-Ji region mainly has a cold temperate climate with dry winters and hot summers, corresponding to the Dwa type in the Köppen-Geiger climate classification [37].

This region provides an ideal study area for this study. First, the tremendous urban expansion has resulted in great changes in the urban form, providing diverse urban samples. Second, the abundant, densely distributed cities with uniform climatic conditions make it easy to control for climate factor variables and biomes, thereby eliminating the effects of regional natural factors on the UHI effect [33]. Third, the temperate climate of this region, which features distinct summers and winters, makes the seasonal variation in the UFI-UHI relationship easier to observe and has considerable reference value for other cities globally. We considered connected urban areas across different administrative units as single research units and identified 150 urban areas for study.

2.2. Data Sources, Metrics, and Processing

2.2.1. Urban Heat Island Intensity

We used monthly composite land surface temperature (LST) remote sensing products at 1-km resolution derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) (Terra) dataset to calculate the surface UHI intensity (SUHII). The SUHII of each city was calculated based on the differences between the LSTs in urban and rural areas [38] using Equation (2):

SUHII = LST_urban − LST_rural

(2)

where SUHII represents the surface UHI intensity and LST_urban and LST_rural represent the LSTs in urban and rural areas, respectively. The urban extents were extracted from the urban category of the land use data with a 30 m resolution acquired from the Resource and Environment Data Cloud Platform, CAS (http://www.resdc.cn/DataList.aspx). The extents of the rural areas were identified based on buffer zones with 10 km radii around these urban extents (shown in Figure 2d–g) [39,40]. Based on the LST data in the years 2000, 2005, 2010, and 2015, we obtained the SUHII datasets of summer day (SUHII_SD), summer night (SUHII_SN), winter day (SUHII_WD), winter night (SUHII_WD), and the annual SUHII (SUHII_YEAR).

2.2.2. Urban Form Measures

We use four categories of indicators to study their relationship with SUHII: urban shape geometric indicators (GEO), urban size and density indicators (SIZ), urban vegetation indicators (VEG), and atmospheric indicators (ATM) (Appendix A

Table A1. Calculation methods and descriptions of the urban form indicators used in this study.

UFI	Equation	Parameter	Description
Urban Area Size (AREA)	$\mathrm{AREA}=\log \left({\mathrm{a}}_{\mathrm{i}}\right)$	ai = area (km2) of urban grid i.	AREA indicates the area of urban expansion.
Urban Population Density (POP)	$\mathrm{POP}=\log \left(\frac{{\mathrm{p}}_{\mathrm{i}}}{{\mathrm{a}}_{\mathrm{i}}}\right)$	pi = population size of urban path i, ai = area size of urban patch i	POP represents the degree of population density.
Nighttime Light Intensity (NLI)	${\mathrm{NLI}}_{\mathrm{ref}}{=\mathrm{aNLI}}_{\mathrm{tar}}{+\mathrm{bNLI}}_{\mathrm{tar}}^2+\mathrm{c}$	NLIref = the reference image NLItar = the target image a,b,c = the coefficients estimated using the ordinary least squares regression technique.	NLI shows the lights generated from electricity. Areas of high economic prosperity and population are generally the areas that are well illuminated.
Fractal Dimension Index (FRAC)	$\mathrm{FRAC}=\frac{2\ln ({.25\mathrm{p}}_{\mathrm{ij}})}{{\mathrm{lna}}_{\mathrm{ij}}}$	pij = perimeter (m) of patch ij. aij = area (m2) of patch ij.	FRAC is used to quantify the degree of complexity of the planar shapes.
Contiguity Index (CONTIG)	$\mathrm{CONTIG}=\frac{[\frac{{\displaystyle \sum _{\mathrm{p}=1}^{\mathrm{n}}{\mathrm{c}}_{\mathrm{ijp}}}}{{\mathrm{a}}_{\mathrm{ij}}^{*}}]-1}{\mathrm{s}-1}$	cijr = contiguity value of pixel p in zone ij. aij* = area of zone ij is represented by the number of grids. s = sum of the values in a 3-by−3 grid template.	CONTIG assesses the spatial connectedness or contiguity of cells within a grid-cell patch. It equals 0 for a one-pixel patch and increases to a limit of 1 as patch contiguity or connectedness increases.
Elongation index (ELONG)	$\mathrm{ELONG}=1-[\frac{{\mathrm{a}}_{\mathrm{ij}}}{{\mathrm{a}}_{\mathrm{ij}}{}^s}]$	aij = area (m2) of grid ij. aijs = area (m2) of smallest circumscribing circle around grid ij.	ELONG is used to distinguish patches that are linear (narrow) or elongated.
Normalized Difference Vegetation Index (NDVI)	$\mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}$	NIR = a subset of the infrared band of the electromagnetic spectrum. Red = the reflectivity corresponding to the red channel.	NDVI is a normalized transform of the NIR to red reflectance ratio, and it is commonly designed to standardize Vegetation Indices values to between −1 and +1.
Tree Cover (TREE)	$\mathrm{TREE}=\frac{{\mathrm{a}}_{\mathrm{tree}}}{{\mathrm{a}}_{\mathrm{sum}}}$	atree = area (m2) of urban forest coverage. asum = area (m2) of the whole urban region.	TREE describes the percent of urban forest coverage.

). Shape continuity (CONTIG) is used to assess the spatial connectedness of urban patches [41]. Fractal dimension (FRAC) is applied to quantify the degree of planar complexity of the urban boundaries [16]. GEO indicators were derived from land-use data. Elongation (ELONG) is represented by the related circumscribing circle index and measures the overall urban patch elongation [42]. Urban area size (AREA) reflects the sizes of urban patches, which are closely related to the area of impervious surfaces. Population density (POP) is calculated with grid data at 1-km resolution obtained from the Center for International Earth Science Information Network at Columbia University (http://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-rev10). These data reflect the population density of the city, and POP has a direct relationship with the purpose of the buildings. The night light index (NLI) was calculated based on NOAA’s nighttime light index data with an approximate 1-km resolution (http://ngdc.noaa.gov/eog/dmsp/) from the Defense Meteorological Satellite Program Operational Linescan System (DMSP-OLS). The NLI can serve as the urban development index (UDI) [43,44]. We calibrated the original NLI data using the ridgeline sampling regression method to obtain a consistent NLI time series [45]. Due to the lack of NLI data in 2015, we used 2013 data for 2015. We used log transforms for urban AREA and POP values to present data closer to a normal distribution to better establish models. The normalized difference vegetation index (NDVI) was derived from MODIS datasets (MOD13Q1) with 500 m resolution and used to estimate the coverage of all types of green plants within the urban areas. Tree cover (TREE) was obtained from MODIS datasets (MOD44B) with 250 m spatial resolution to measure tree cover. We calculated the difference between the urban and rural areas (annotated as dNDVI and dTREE, respectively) to better model the relationships [15]. A larger dTREE value indicates larger forest areas, which provide more shade and generate more transpiration during the day than other, shorter plants. The descriptions of all the UFIs are shown in Appendix A Table A1.

2.2.3. Atmospheric Conditions

We considered several atmospheric factors (ATM) as control variables to evaluate their influences on SUHIIs [46]. The ATM factors are annual and seasonal wind speed (WIN), relative air humidity (RHU), amount of precipitation (PRE), and aerosol optical depth (AOD). The first three indicators are derived from daily meteorological monitoring data acquired from the China National Meteorological Center (http://cdc.cma.gov.cn/) and are interpolated to a raster format using the cokriging interpolation method. The AOD data are derived from the MODIS dataset (MOD04C6) and have a 10-km spatial resolution (https://modis.gsfc.nasa.gov/data/dataprod/mod04.php). All the raster layers with spatial resolutions larger than 1 km were resampled to 1 km for zonal statistics.

2.3. Statistical Models

2.3.1. Multivariate Ridge Regression (M-RR)

The widely used ordinary least squares (OLS) method can fail to assign proper weights or signs to the individual explanatory variables used as predictors in many cases [47,48]. As an alternative, multivariate autoregressive (M-AR) models can be applied. However, such discrete processes (whose variables are either retained or discarded) often exhibit a high variance and result in the incomparability of coefficients of the selected and unselected variables in the different models [49]. Here, we use a multivariate ridge regression model (M-RR) to avoid the aforementioned distortions and to estimate the coefficients of each independent variable for comparison [50,51]. Assuming that SUHII_i is the SUHII of the ith city (i = 1, 2, 3, …, N), UFI_ij is the jth UFI of the ith urban area, β₀ is the bias, and β_j is the coefficient of the jth UFI, the description of RR is as follows:

$${\stackrel{\wedge }{\beta }}^{\mathit{ridge}}=\underset{\beta }{argmin}\left\{{\displaystyle \sum }_{i=1}^N(SUHI{I}_i-{\beta }_0-{\displaystyle \sum }_{j=1}^{p}UF{I}_{ij}{\beta }_j{)}^2+\lambda {\displaystyle \sum }_{j=1}^p{\beta }_j^2\right\}$$

(3)

where λ is a complexity parameter controlling the amount of shrinkage. By imposing λ as a size constraint on the coefficients, we impose very large positive coefficients on one variable, which are canceled by the similarly large negative coefficients on its correlated partner. We adopted the generalized cross-validation (GCV) method to obtain a more reliable estimate of λ [52,53,54].

2.3.2. Other Models

As a comparison, all-subset M-AR models and panel data models were used [55,56]. The adjusted R² value was taken as the indicator for judging the fitness of an M-AR model. Due to the limited number of time periods, univariate panel models were constructed [20]. The Hausman test was used to decide whether the univariate fixed effects (U-FE) or univariate random effects (U-RE) model should be used [57]. We defined the relative contribution as the proportionate contribution each predictor makes to the R² value considering both the unique contribution of each predictor by itself and its incremental contribution when combined with the other predictors [58]. Figure 3 shows the overall framework of the data processing and statistical analysis performed in this study.

3. Results

3.1. Seasonal and Diurnal Variations in the SUHII Distribution

The spatial distribution of SUHIIs showed significant seasonal and diurnal variations. In summer, the SUHII in the southern Jing-Jin-Ji region was relatively low during the day but increased during the night (Figure 4b,c). In winter, the nighttime SUHIIs were generally more intense than the daytime SUHIIs (Figure 4d,e). The SUHII data distribution showed notable differences among the four time periods (Figure 4f). The summer SUHII was generally higher than that in winter. The SUHII was the lowest in the winter daytime, probably due to the reduction in the incoming solar shortwave radiation-caused heating in this period [59].

3.2. Correlations between the UHIs and UFIs

Figure 5 shows the strong (ρ ≥ 0.2) and weak or no (ρ < 0.2) correlations between UFIs and the four types of SUHIIs. The correlation varied greatly for the four periods. For daytime SUHII, all indicators showed significant seasonal variations (i.e., the difference between the summer ρ and winter ρ was greater than 0.2), except for CONTIG. For the nighttime SUHII, only dNDVI and dTREE showed significant seasonal variations. In both seasons, all eight urban form factors exhibited significant diurnal differences.

3.3. Estimation of the UFI Effects

We constructed three univariate regression models and six multivariate regression models for comparison: univariate OLS regression (U-OLS), univariate fixed effect panel model (U-FE), univariate random effect panel model (U-RE), multivariate all-subset OLS regression (M-AR), three-variable OLS regression (T-OLS), multivariate fixed effect panel model (M-FE), multivariate random effect panel model (U-RE) and M-RR model. The M-AR resulted in the best model (M-BAR) (Appendix A Figure A1). The best λ values of the RR models were obtained by the GCV method shown in Figure 6a and were used to construct the best multivariate ridge regression model (M-BRR). The results of the Hausman test and univariate models are summarized in Appendix A

Table A2. Results of regression models and the Hausman test.

Time	UFI	Univariate OLS Linear Model		Fixed Effects Model		Random Effects Model		Hausman Test
		Coefficient	Signif.	Coefficient	Signif.	Coefficient	Signif.	Chisq	Signif.
SD	CONTIG	4.98		13.32	*	10.92	*	0.98	*
	FRAC	14.85	***	5.84	***	7.91	***	14.15	***
	ELONG	4.26	***	−0.14		1.13	*	28.71	***
	NLI	4.11 × 10−2	***	3.73 × 10−4	***	3.83 × 10−4	***	0.28
	AREA	0.88	***	1.26	***	1.13	***	3.01	*
	POP	0.38	***	3.86	***	0.61	***	58.31	***
	dNDVI	−7.42	***	−9.67 × 10−4	***	−8.57 × 10−4	***	5.98	*
	dTREE	−13.30	***	−0.15	***	−0.14	***	0.71	*
SN	CONTIG	24.75	***	12.13	***	18.73	***	8.61	***
	FRAC	2.73	***	0.80		1.38		3.09	*
	ELONG	−0.90	***	−0.45		−0.61	*	1.35	*
	NLI	2.23 × 10−2	***	1.87 × 10−4	***	2.02 × 10−4	***	1.67	*
	AREA	0.63	***	0.68	***	0.65	***	0.24
	POP	0.37	***	2.14	***	0.43	***	69.17	***
	dNDVI	8.73 × 10−1	***	−1.89 × 10−4	***	−2.48 × 10−5		35.04	*
	dTREE	1.83	***	−0.07	***	−0.01	*	79.23	***
WD	CONTIG	1.49		5.02		1.98		0.47	*
	FRAC	1.56		−2.60		0.76		8.37	**
	ELONG	0.51		−1.59	**	0.19		15.28	***
	NLI	9.12 × 10−3	**	3.48 × 10−5		8.33 × 10−5	**	1.61	*
	AREA	0.16	*	−0.08		0.14		2.38	***
	POP	0.06		0.67		0.06		1.48	*
	dNDVI	1.68	**	4.55 × 10−4	***	2.29 × 10−4	***	11.36	***
	dTREE	−3.44	***	0.02		−0.03	***	6.82	**
WN	CONTIG	32.90	***	0.88		9.78	***	130.28	***
	FRAC	2.38	*	2.99	***	2.91	***	0.15
	ELONG	−1.93	***	0.65	**	0.21		38.29	***
	NLI	2.53 × 10−2	***	1.98 × 10−4	***	2.06 × 10−4	***	2.04	*
	AREA	0.79	***	0.41	***	0.52	***	12.91	***
	POP	0.46	***	0.71	***	0.48	***	1.38	*
	dNDVI	−1.73	**	−2.44 × 10−4	***	−2.35 × 10−4	***	0.63	*
	dTREE	7.81	***	−0.04	***	0.02	**	154.24	***

Note: * is significant at the 0.05 level (2-tailed); ** is significant at the 0.01 level (2-tailed); and *** is significant at the 0.001 level (2-tailed).

. The estimated slope coefficients for all models are summarized in

Table 1. Regression coefficients and fitness of the models.

Time	Model	Slope Coefficient								Performance
		CONTIG	FRAC	ELONG	NLI	AREA	POP	dNDVI	dTREE	R2	R2 (CV)
SD	M-BRR	16.32	4.42	1.11	0.02	0.06	0.42	−6.44	−3.88	0.5014	0.4259
	M-BAR	22.56	5.53	1.18			0.40	−6.54	−4.01	0.5014
	M-OLS				11.42		2.00	−8.19		0.5015	0.4256
	T-OLS		4.87			0.81		−7.33		0.4503
	U-OLS		14.85	4.26	0.41	0.88	0.38	−7.42	−13.30	0.1327
	U-FE	13.32	5.84		0.00	1.26	3.86	−9.67	−15.20
	U-RE	10.92	7.91	1.13	0.38	1.13	0.61	−8.57	−14.06
	M-FE			−1.11	0.12		1.56	−1.58	−0.07	0.2445
	M-RE				0.12	0.31	0.28	−2.17	−12.20	0.3037
SN	M-BRR	3.98	0.38	−0.14	0.04	0.10	0.06	0.13	0.25	0.3545	0.2372
	M-BAR	11.00			0.10	0.11	0.15	0.45		0.3735
	M-OLS			−0.48	0.07	0.33	1.10		−3.54	0.3749	0.2185
	T-OLS					0.66		0.70		0.2804
	U-OLS	24.75	2.73	−0.90	0.22	0.63	0.37	0.87	1.83	0.1370
	U-FE	12.13			0.00	0.68	2.14	−1.89	−7.21
	U-RE	18.73	1.38	−0.61	0.20	0.65	0.43		−1.40
	M-FE			−0.48	0.06	0.28	1.04	−1.19	−2.67	0.2956
	M-RE	6.27		−0.47	0.09	0.27	0.14	−1.98		0.3245
WD	M-BRR	0.02	−0.20	−0.01	0.05	0.05	0.07	1.96	−3.20	0.0830	0.0646
	M-BAR				0.04		0.12	2.58	−4.07	0.0859
	M-OLS			−1.28	0.12		1.69	5.98		0.0868	0.0568
	T-OLS					0.15		1.92		0.0222
	U-OLS				0.09	0.16		1.68	−3.44	0.0116
	U-FE			−1.59				4.55
	U-RE				0.08			−2.29	−3.27
	M-FE			−1.28	0.12	1.69	0.00			0.0914
	M-RE							2.96	−4.17	0.0754
WN	M-BRR	0.95	−0.55	−0.39	0.18	0.39	0.01	−2.02	7.39	0.5690	0.5075
	M-BAR	5.16			0.20	0.30		−2.21	8.00	0.5702
	M-OLS				18.17		−0.58	−1.00		0.5710	0.4985
	T-OLS		−4.06			0.89		−1.46		0.2576
	U-OLS	32.90	2.38	−1.93	0.25	0.79	0.46	−1.73	7.81	0.1580
	U-FE		2.99	0.65	0.00	0.41	0.71	−2.44	−4.27
	U-RE	9.78	2.91		0.21	0.52	0.48	−2.35	2.03
	M-FE				0.18		−0.58	−1.00		0.2392
	M-RE	4.60			0.19	0.14		−1.18	5.42	0.3336

Note: M-BRR = best multivariate ridge regression model, M-BAR = best all-subset OLS regression model, M-OLS = multivariate OLS regression model, T-OLS = three-variable OLS regression model, U-FE = univariate fixed effect panel model, U-RE = univariate random effect panel model, M-FE = multivariate fixed effect panel model, and M-RE = multivariate random effect panel model. The blank entries in the table indicate that the variables were not used, or the tests were not performed. The green crosses indicate that the variable did not pass the significance test at the 0.1 level. R-Square (CV) = mean R-square value of the three-fold cross-validations. In addition, the NLI coefficient has been multiplied by 10 for demonstration purposes, and the R² of U-OLS is the average of eight values.

. In the table, we ignored those slope values with a confidence less than 90%, except for M-RR, whose slope coefficients had already been degraded by the penalty term. Our results showed poor fitness of all models during WD but good fitness of the models at other times. The RR coefficients were rescaled by dividing by the maximum absolute value. Thus, we could compare the difference between the effects of the same UFI on SUHIIs across different time periods, as shown in Figure 6b.

In general, the results for all univariate models, including panel models, indicated poor fitness, with the slope values fluctuating sharply. In addition, the Hausman test results for each variable were contradictory for different models. These results indicate that the panel models with few explanatory variables are inferior to the M-RR model in solving the problem of missing variables. Among the multivariate models, although M-AR had a good explanatory power, it is difficult to compare the slope values of the different best models for each UFI due to the discarded variables. Considering robustness confirmed by the three-fold cross-validation results (R-squared (CV)) (Table 1), we can conclude M-RR as the best model to estimate the effects of all UFIs.

There were notable diurnal and seasonal variations among the slope coefficient signs of the three specific GEO factors. The signs of CONTIG were always positive in the four models, whereas the signs of FRAC were opposite in winter and summer, and those of ELONG were opposite in the daytime and nighttime. The absolute slope values of most GEO factors were larger in summer than in winter. In particular, an increase in CONTIG resulted in a four-fold stronger SUHII increment in the summer daytime than in the summer nighttime. The same FRAC increment enhanced the SUHII intensity by approximately 10-fold in the summer daytime than in the summer nighttime, while FRAC weakened the SUHII effect in winter with an influence magnitude similar to that in the summer nighttime. An increased ELONG would result in an increased daytime SUHII and decreased nighttime SUHII. The magnitude of its influence was greatest during the summer daytime.

The regression coefficient signs of the SIZ factors were all positive. However, the slope values for NLI and AREA were higher in winter than in summer and higher at night than during the day. Furthermore, the variance pattern of POP was similar to that of the GEO factors, with higher slope values in summer than in winter and higher values in the daytime than at night. In summer, dNDVI and dTREE presented significant daytime cooling effects, which were tens of times stronger than the night temperature preservation effect. However, the situation became more complex in winter, as the slope values of dNDVI were positive in the daytime and negative at night, whereas the pattern of dTREE was the opposite. Moreover, the absolute values of the slope of dTREE in the nighttime were much larger than those in the daytime, while the two slope values of NDVI were almost equal.

3.4. Contributions of Three Main Urban Form Categories

We initially calculated the relative contribution rates of each urban form index and meteorological factor (Appendix A

Table A3. Relative contribution rate.

	CONTIG	FRAC	ELONG	NLI	AREA	POP	dNDVI	dTREE	AOD	PRE	RHU	WIN
SUHII_YEAR	12.28%	5.75%	1.76%	21.70%	15.94%	15.69%	5.63%	7.45%	2.10%	5.05%	2.85%	3.80%
SUHII_SD	3.35%	5.16%	5.67%	5.36%	3.80%	6.81%	31.61%	23.51%	2.45%	4.05%	7.38%	0.85%
SUHII_WD	1.66%	1.44%	1.85%	4.78%	1.95%	3.52%	12.13%	36.07%	7.22%	2.16%	23.29%	3.93%
SUHII_SN	20.66%	3.22%	1.97%	21.66%	15.32%	19.56%	1.56%	1.84%	2.77%	3.89%	6.22%	1.32%
SUHII_WN	10.37%	1.98%	3.66%	14.85%	12.28%	8.24%	2.15%	31.36%	13.82%	0.85%	0.36%	0.09%

). By summing the contributions of all indicators in each dimension, we obtained the total contributions of the three dimensions of the urban form and atmospheric factors (Figure 7a). The results showed that VEG dominated the daytime SUHII, contributing 55.12% in the summer daytime and 48.20% in the winter daytime. SIZ dominated the nighttime SUHII, contributing 56.54% and 35.37% in the summer daytime and in the winter daytime, respectively. Atmospheric factors, including air pollution, played a large role in winter, similar to that of vegetation. At the annual scale, the contribution rate of SIZ was the highest, followed by the GEO, and finally by the ATM and VEG.

By summing the contribution rates of each of the four periods in a year, we identified the trend of the total contribution of the four urban form categories over decades (Figure 7b). In general, there is no obvious seasonal or diurnal difference in the trend of contribution rates for various urban forms. The contributions of GEO and ATM increased continually, and that of VEG decreased continually.

4. Discussion

4.1. Comparison Overview of the Correlation, Slope Coefficient, and Relative Contribution Analyses

Previous studies on UFIs and SUHIIs have mostly focused on the correlation coefficient or goodness of the regression fit rather than the regression slope coefficient. Although both the correlation coefficient and goodness of the regression fit can indicate the degree of correlation, the regression slope coefficient provides a quantitative tool for comparing the UFIs’ impacts across multiple time periods. Taking the SIZ indicators as an example, we find that their correlations with SUHIIs are higher at night, but their slope coefficients are higher during the day. These results are consistent with those of many previous studies [10,33,60,61]. However, in these studies, the estimation of slope coefficients did not contribute to the conclusions. For example, a study of 65 select cities in North America found that the correlation between the population size and SUHII was higher during the night than during the day (with R² values of 0.54 and 0.20, respectively) but ignored the fact that the slope coefficient was much higher during the daytime than during the nighttime (1.60 and 0.64, respectively) [10]. Similar results can also be found in studies by Tran et al. [60], Tan and Li [33], and Li et al. [61]. The results of all of the above studies implied that although there is a slightly stronger correlation between population density and SUHII at night, the rise in population density has a greater impact on SUHII intensity during the day.

Currently, no unanimous conclusion has emerged on the seasonal or diurnal variations in the relative contributions of urban form factors [62,63]. The lack of quantitative attribution of the various contributions to SUHIIs has been considered a serious barrier to UHI mitigation due to the difficulty of obtaining the most efficient scheme with limited manpower and material resources [10,39,64,65]. Through relative contribution rate estimation, we found that there were significant seasonal and diurnal variations in terms of the dominant urban form factor. Moreover, their contributions throughout a year have changed regularly over the years. The results detailed in Section 3.4 are partially similar to a study of 193,090 cities globally, whose rank in terms of relative importance can be listed as EVI, AREA, NLI, and POP during the daytime and AREA, NLI, EVI, and POP during the nighttime [39].

4.2. Reasons for the Seasonal and Diurnal Variations in the Urban Form Effects

FRAC shows the opposite seasonal pattern, which is supported by existing studies. For instance, Tu et al. [66] proved that there was a negative correlation between FRAC and SUHII in October. However, Zhou et al. [16] used summer LST to prove that FRAC and UHI intensity are positively correlated. This seasonal difference is due to the multiple effects of FRAC. On the one hand, a high FRAC value demonstrates that the city has a large circumference–area ratio, which means that the urban areas and their peripheral areas exhibit a high degree of contact and that fresh, cool air can easily enter from the surrounding rural areas. On the other hand, a high FRAC value is also associated with a greater traffic burden and high heat and greenhouse gas emissions [18,20,67], which enhance the UHI effect. In addition, cities with high FRAC values feature intertwined built-up and non-built-up areas in the fringe zone, which can increase the surface aerodynamic roughness. Notably, this condition is also context-dependent and is closely related to the specific land-use characteristics. Due to the seasonal and diurnal variations in energy absorption and release, the two opposite effects together produce various overall effects. The impact of ELONG on SUHII is opposite between day and night, and thus, the correlation between ELONG and SUHII_YEAR is not significant.

The positive and negative effects of SIZ factors on SUHII are consistent among different seasons and times of day but their regression slope coefficients have significant seasonal and diurnal variations. This finding is consistent with those of previous studies showing that the slope coefficients of NLI and AREA are both greater in winter than in summer [10,33,68]. Variation patterns of NLI were also similar to the findings in a study of the Yangtze River Delta urban agglomeration [68]. The reason may be that larger and more economically prosperous cities emit more heat in winter and at night. The additional urban energy sources during winters and nights can significantly affect the UHI intensity. The effect of AREA is greater at night, which may be due to the increased heat capacity caused by the urban expansion and the resulting increase in nighttime heat radiation. The slope coefficient of POP in summer was much higher than that in winter and higher in the daytime than at night, which was supported by other studies [15,68]. This situation may partly exist because the buildings are denser in more densely populated areas, decreasing ventilation. Changes in these built environment characteristics were found to be more influential in hotter seasons [63].

The positive or negative impacts of dTREE and dNDVI were the same in summer but opposite in winter. This heterogeneity in the vegetation effect has been confirmed by previous studies and is considered to be highly context-dependent, particularly with respect to geographical conditions and climatic backgrounds [10]. In our context, the NDVI is a measure of the coverage of all green plants, while the TREE index measures the continuous forest coverage. Compared with general plants, forest trees exhibit greater transpiration, shading, and wind resistance. In summer, the transpiration of vegetation is strongest in the daytime, resulting in a strong cooling effect. However, during summer nights, transpiration decreases, and the wind-blocking effect of plants becomes more prominent. Therefore, in summer, dNDVI and dTREE were negatively correlated with SUHII during the day and positively correlated with SUHII at night. In winter, the transpiration of deciduous plants in the region is weak but the plants help absorb heat by reducing the surface albedo, thereby enhancing the daytime SUHII [15]. However, at night, plants can generate mechanical turbulence and enhance the convective heat loss, thereby helping to mitigate SUHII [10,38]. A larger dTREE value indicates larger forest areas, which provide more shade and generate more transpiration during the day than other, shorter plants. Therefore, a high dTREE value provides a slight cooling effect in winter daytime. However, at night, the wind resistance caused by the presence of densely distributed trees becomes prominent, thereby increasing the UHI intensity.

4.3. Urban Planning Implications

Although there is a consensus that compact cities worldwide are preferable in the context of sustainable development, debates remain about what urban form is beneficial to the thermal environment at the city scale [13,17]. The seasonal and diurnal variations in the relationship between urban form and the UHI effect in this paper provide the following insights into scientific urban planning. First, to effectively mitigate the impact of summer heatwaves, UFIs that play a greater role in summer should receive the most attention. Overly continuous and overly elongated cities with highly complex boundaries or high population densities should be avoided in urban planning. Second, because of the cold winters in this area, an urban form that can help alleviate the UHI effect in summer and enhance the UHI effect in winter would help conserve energy [32]. For example, afforestation in cities, rather than simply increasing the green space coverage through grassland or short vegetation, is one such win-win solution. Our results showed that the all-day average SUHII will increase 0.42 °C in winter for each 0.1 increase in dTREE, while the opposite effects of a dNDVI increase would almost equally offset this increase in this region. Third, although it seems to be necessary to control the scale of cities in high-heat-risk urban areas because of its consistent effects diurnally and seasonally, additional attention should be paid to optimization of the urban geometric form. This geometric optimization need not be at the expense of urban development, and its relative contribution increased over the 15-year study period.

4.4. Limitations and Future Studies

The first limitation is that most of our research samples are small- and medium-sized cities. Future research can expand the observation scope to larger cities under different climatic backgrounds around the world to provide more comprehensive insights. Second, as an empirical study, more detailed explanations of the mechanism of our urban-scale statistical results are beyond the scope of this paper. More convincing interpretations can be obtained in combination with urban climate models. Third, our models did not explain the winter daytime SUHIIs well. However, the RR models with penalty terms have minimized the bias. In addition, more measurements of the urban form, such as building height, building density, water area, and the mix of land use, etc. would be useful. Researchers can supplement these measurements of the urban form based on geographic information big data. At the city scale, these indicators would help to provide more information about the city’s social and economic activities and provide a better perspective on the impact of ventilation and other aspects.

5. Conclusions

Our results show that urban form has an important influence on the UHI effect. However, there are significant seasonal and diurnal differences in the impacts of all urban form indicators. First, all the urban size factors enhance the UHI effect, but urban geometry and vegetation factors can have both positive and negative impacts. The relationships between urban shape complexity, elongation, NDVI, and forest coverage and the UHI effect show various patterns depending on the season and time of day. Second, there are significant seasonal and diurnal variations in the magnitude of the influence of urban form, as measured by the regression slope coefficients. The influences of urban continuity, complexity, elongation, and population density in summer are significantly higher than those in winter, and those during summer daytime are higher than those during summer nighttime. The impacts of NLI and AREA on SUHII are greater in winter than in summer and at night than during the day. The effect of dNDVI during summer daytime is significantly higher than that in other periods. Third, in terms of relative contribution, urban vegetation factors dominate during the day, and size factors dominate during the night. Urban geometry can be a subdominant factor in summer, but its contribution is much smaller in winter. The relative contributions of the geometric factors increased consistently over the 15-year study period, those of the scale factors decreased over the last five years of the study period, and those of the vegetation factors decreased consistently.

To improve the urban thermal environment, three major urban planning proposals have been proposed. First, it is necessary to avoid overly continuous, disorderly sprawl and overly elongated cities. This measure is important to mitigating the impact of summer heatwaves. Second, afforestation in cities, rather than simply expanding the coverage of short vegetation such as grassland, is of great value to alleviating the cooling or heating energy burden across all seasons. Third, in addition to the size of the city, more attention should be paid to the optimization of geometric form in the future because the relative contribution of geometric form has increased over time and such optimization need not be detrimental to urban economic development.