An Assessment of Urban Surface Energy Fluxes Using a Sub-Pixel Remote Sensing Analysis: A Case Study in Suzhou, China

Abstract

Urban surface energy fluxes are closely associated with land-cover types (LCTs) and critical biophysical compositions. This study aims to assess the contribution of LCTs, vegetation fractional coverage (VFC) and percentage of impervious surface area (ISA%) to urban surface energy fluxes using remote sensing. An advanced urban surface energy flux algorithm was used to combine satellite imagery and meteorological station data to investigate the thermal environments in the city of Suzhou, China. The land cover abundances retrieved by multiple endmember spectral unmixing analysis (MESMA) were used to retrieve the per-pixel sensible heat flux (H) and latent heat flux (LE). The resultant heat fluxes were assessed using evaporation pan data collected from meteorological stations and ratios of the heat fluxes to the net radiation (Rn). Furthermore, spatial patterns of urban heat energy were investigated using an integrated analysis among land surface temperature (LST), heat fluxes, LCTs, VFC and ISA%. The high values of H and LST were found over the urbanized areas, which also had low values of LE. Conversely, the vegetated area was characterized with high LEs, as well as low LSTs and Hs. Moreover, a statistically-significant correlation (p < 0.05; R² = 0.88) was observed between LE and VFC at the zonal level, and a statistically-significant correlation (p < 0.05; R² = 0.90) was exhibited between H and ISA%. It is concluded that VFC, ISA% and LCTs are promising for delineating urban heat fluxes. Overall, this study indicates that remote sensing techniques can be used to quantify urban thermal environments.

1. Introduction

The urban heat island (UHI) effect generally results in increased local temperatures in urban areas compared to the surrounding rural areas. This effect is due primarily to changes resulting from artificial land surfaces and ongoing human activities. Surface temperatures and urban energy fluxes have been studied by many researchers to investigate the UHI. Moreover, urban thermal environments are closely associated with the land cover types (LCTs), especially impervious surfaces and vegetation cover [1,2,3,4,5]. Impervious surfaces and vegetation cover are two critical components that modify urban surface energy fluxes and eventually lead to differences in urban heat energy between the city and the countryside. Through the characteristics of impervious materials and vegetation cover, urban surfaces can influence energy exchange processes, such as solar radiation absorption and thermal emission, as well as latent and sensible heat fluxes. Since [6] found that a higher level of latent heat flux (LE) was correlated with vegetated areas and that sensible heat flux (H) was more favored by urban impervious areas, studies have been conducted to explore the contributions of impervious surfaces and vegetation to surface energy exchange [3,7,8]. Furthermore, it has been reported that combining statistical approaches and energy balance models can achieve an enhanced characterization of the patterns of surface landscape elements, thereby facilitating future applications of urban thermal environment analyses. On the other hand, although the concepts of landscape elements and urban thermal energy fluxes have been known for some time, a quantitative interpretation of their roles and interactions across complex urban environments remains elusive. A better understanding of the effects of LCTs and two biophysical parameters (vegetation fractional coverage (VFC) and percentage of impervious surface area (ISA%)) on LSTs and energy fluxes is required in urban studies.

The ground measurement-based method and the remote sensing method are two techniques that are broadly applied to explore urban surface energy fluxes. Ground measurement-based studies can delineate detailed spatial-temporal variations in urban thermal patterns [9,10,11]. However, as an independent technique constrained by physical and economic factors, this method has been employed only over small local-scale areas. On the other hand, a diversity of satellite remote sensing data has been extensively used to explore urban thermal patterns over a wide area [12]. Specifically, for analyzing the energy budgets of urban regions, remote sensing data have been coupled with numerical models to enhance modeling accuracy. For example, [13,14,15] attempted to study urban heat energy fluxes based on a heat-balance model using medium spatial resolution remote sensing and ground meteorological data. Xu et al. [16] used hyperspectral imagery from the Operative Modular Imaging Spectrometer, a survey map and meteorological data to map spatial variations in turbulent sensible heat flux in Shanghai, China. However, concerning the resolution of the most used operational satellite sensors (i.e., LANDSAT and ASTER), the limitations in the spatial dimension for characterizing urban heat energy are evident. Urban environments are generally accompanied by a mixture of impervious surfaces and vegetation cover within a satellite sensor pixel. It should be noted that most urban studies have focused on describing surface heat fluxes based on traditional hard classification schemes [14,17,18]. The utility of biophysical composition parameters (i.e., sub-pixel-level VFC and ISA%) for the assessment of urban energy budgets has not received the attention it deserves, in particular at local scales, in spite of the existence of a few interesting studies on the thermal environment (e.g., [19,20,21,22]).

The vegetation-impervious surface-soil (VIS) model, which is among those models that are used to study the biophysical composition of urban environments with remote sensing spectral characteristics [23], provides an alternative approach for solving the pixel mixture problem, especially for complex urban landscapes. This conceptual model assumes that surface pixels can be expressed by a combination of three components: vegetation, impervious surfaces and soils [24,25]. The V-I-S model has been widely implemented using the approach of spectral mixture analysis (SMA) in urban environment applications. Furthermore, biophysical compositions obtained using SMA have reportedly achieved success in assessing the contributions of landscape elements to urban thermal environments [1,22,26]. In spite of this, the related research that has focused on using the SMA approach to study urban surface energy fluxes, especially regarding improving estimated urban heat fluxes and analyzing their interactions with urban biophysical compositions, seems insufficient.

This study has attempted to investigate spatial patterns in the LSTs, surface energy fluxes and the interplay with ISA% and VFC in the city of Suzhou, China. In particular, multiple endmember spectral mixture analysis (MESMA), which was presented in [27] and has been extensively applied in a variety of fields [28,29], was used to retrieve accurate VFC and ISA% to effectively estimate per-pixel latent heat flux and sensible heat flux. Overall, the specific objectives of this research are (1) to describe the spatial patterns in urban heat energy (LST and heat fluxes) and (2) quantitatively examine the relationships among LST, heat fluxes and urban biophysical compositions (VFC and ISA%).

2. Materials and Method

2.1. Study Area and Materials

An area of approximately 1200 km² that covers the city of Suzhou, China, was selected for the present research. Suzhou is located in eastern China (Figure 1) and has a population of over 10 million. The center of the present study area is the urban areas of Suzhou, and the remaining area is mostly used for agriculture and water storage. Previous studies [30,31] have reported that Suzhou has experienced significant urban expansion in recent decades. Examining the thermal environmental impacts of urban sprawl in Suzhou is beneficial for understanding and planning its development and is especially conducive to mitigating the adverse thermal effects of the layout of the urban surface landscape.

Remotely-sensed satellite data and ground meteorological data were collected for this research. LANDSAT-5 TM (L1G-level product, path/row 119/038) data for 24 May 2010 were acquired. The L1G data were geometrically corrected only to remove systematic geometrical errors. The image provides six visible and near-infrared bands and a thermal infrared band, the spatial resolutions of which are 30 m and 120 m, respectively. The LANDSAT-5 TM image was rectified to the Universal Transverse Mercator coordinate system based on a high resolution Google Earth™ image, and a root mean square error of less than 1 pixel was obtained for the rectification. Atmospheric corrections were applied to the visible and near-infrared bands of the image using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) module [32]. The ground meteorological data were provided by the Central Weather Bureau of Suzhou. Table 1 shows the locations of the meteorological stations. The data used include the sunshine duration (SSD), wind speed (WS), relative humidity (RH) and atmospheric temperature (AT) acquired on 24 May 2010.

Figure 1. The left map is Landsat-5 TM true color composition imagery acquired over Suzhou city; the top right is the location of Suzhou in China, and the bottom map is for the meteorological stations, the number of which corresponds to that listed in Table 1.

Figure 1. The left map is Landsat-5 TM true color composition imagery acquired over Suzhou city; the top right is the location of Suzhou in China, and the bottom map is for the meteorological stations, the number of which corresponds to that listed in Table 1.

Table 1. Summary of the meteorological conditions (at a 2-m height) at 1040 local time on 24 May 2010 for Suzhou. AT, atmospheric temperature; WS, wind speed; RH, relative humidity; SSD, sunshine duration.

ID	Meteo-Stations	Latitude (degree)	Longitude (degree)	AT (°C)	WS (ms−1)	RH (%)	SSD (Hour)
1	Suzhou	31.27	120.63	24.7	2.2	46.0	11.2
2	Dongshan	31.07	120.43	22.9	2.0	56.0	12.2
3	Wujiang	31.13	120.62	23.0	4.5	54.7	11.6
4	Suzhou Park	31.29	120.54	25.0	3.1
5	Xukou Marine Firm	31.23	120.47	23.1	3.5
6	Xiangcheng Road	31.37	120.63	24.5	3.0	42.7
7	Kuatang	31.33	120.68	24.5	2.6
8	Experimental Primary School	31.28	120.62	23.9	1.2
9	Sanmadun	31.25	120.78	23.2	4.7
10	Huangdi	31.42	120.6	24.6	3.9

Table 1. Summary of the meteorological conditions (at a 2-m height) at 1040 local time on 24 May 2010 for Suzhou. AT, atmospheric temperature; WS, wind speed; RH, relative humidity; SSD, sunshine duration.

ID	Meteo-Stations	Latitude (degree)	Longitude (degree)	AT (°C)	WS (ms−1)	RH (%)	SSD (Hour)
1	Suzhou	31.27	120.63	24.7	2.2	46.0	11.2
2	Dongshan	31.07	120.43	22.9	2.0	56.0	12.2
3	Wujiang	31.13	120.62	23.0	4.5	54.7	11.6
4	Suzhou Park	31.29	120.54	25.0	3.1
5	Xukou Marine Firm	31.23	120.47	23.1	3.5
6	Xiangcheng Road	31.37	120.63	24.5	3.0	42.7
7	Kuatang	31.33	120.68	24.5	2.6
8	Experimental Primary School	31.28	120.62	23.9	1.2
9	Sanmadun	31.25	120.78	23.2	4.7
10	Huangdi	31.42	120.6	24.6	3.9

2.2. Methodology

To assess the spatial pattern of the urban heat energy, an LST map of Suzhou at 1040 local time (LT) on 24 May 2010 was first retrieved using a generalized single-channel algorithm developed by [33]. The heat flux components (Rn, H, LE and storage heat flux) were then estimated based on an improved two-source energy balance algorithm. In addition, a VFC, ISA% and LCTs map was estimated using the approach of SVM and MESMA. The overall flowchart is given in Figure 2.

Figure 2. The overall flowchart of this research; the dashed boxes from top to bottom represent the section of input data, data processing and output data, respectively.

Figure 2. The overall flowchart of this research; the dashed boxes from top to bottom represent the section of input data, data processing and output data, respectively.

2.2.1. Retrieving the LST

For retrieving LSTs using LANDSAT-5 TM data, the thermal infrared band subsequent to radiometric calibration and atmosphere correction was first converted to at-satellite brightness temperatures $T_{sensor}$and then converted to land surface temperature $T_s$. Therefore, we have:

$$T_s=\gamma [{\epsilon }^{-1}({\psi }_1{L}_{sensor}+{\psi }_2)+{\psi }_3]+\delta$$

(1)

$$\gamma ={\{\frac{c_2{L}_{sensor}}{T_{sensor}^2}[\frac{{\lambda }^4}{c_1}{L}_{sensor}+{\lambda }^{-1}]\}}^{-1}$$

(2)

$$\delta =-\gamma L_{sensor}+T_{sensor}$$

(3)

where $L_{sensor}$ is the at-sensor radiance is$W{m}^{-2}s{r}^{-1}$. $\mathsf{\epsilon }$ is the surface emissivity, $\mathsf{\lambda }$ is the effective wavelength, $C_1=1.19104\times {10}^{8}W\mu m^4{m}^{-2}s{r}^{-1}$ and $C_2=14387.7\mu mK$. ${\phi }_1$, ${\phi }_2$ and ${\phi }_3$ are atmospheric functions, which can be calculated using the following expressions:

$${\psi }_1=0.14714w^2-0.15583w+1.1234,$$

(4)

$${\psi }_2=-1.1836w^2-0.37607w-0.52894$$

(5)

$${\psi }_3=-0.04554w^2+1.8719w-0.39071$$

(6)

where $w$ is the water vapor content retrieved from the local meteorological data.

The surface emissivity $\epsilon$ is a key factor in determining LSTs. The Normalized Difference Vegetation Index (NDVI) thresholds method (NDVI^THM) was used to calculate pixel-level emissivity in this study. NDVI^THM [34,35,36] assumes that the vegetation cover and bare soil have specific emissivities and that they can be used to mix other surface pixels. This method obtains the pixel-level emissivities considering different cases. In the case of NDVI <0.2, the pixel is considered to be bare soil, and a known emissivity value of bare soil is then assigned. In the case of NDVI > 0.5, the pixel is considered to be fully vegetated, and the emissivity is obtained from a known emissivity value of vegetation. In the case of 0.2 < NDVI < 0.5, the pixel is composed of a mixture of vegetation and bare soil according to the fractional cover. In this study, the bare soil emissivity was 0.978, and the vegetation emissivity was 0.985 [35].

2.2.2. Retrieval of Heat Flux Components

The surface heat fluxes were estimated following [37] and are represented by Equation (7), in which the absorbed net radiation balances the outgoing fluxes of the sensible heat, latent heat and ground heat:

$$R_n+A=G+LE+H$$

(7)

where $R_n$ is the net radiation, $A$ is the total anthropogenic heat discharge, $G$ is the ground heat flux, $LE$ is the latent heat flux and $H$ is the sensible heat flux.

The net radiation $R_n$ is the sum of the short-wave and long-wave radiation at the land surface, which is calculated as:

$$R_n=\left(1-\alpha \right)K\downarrow +\epsilon {\epsilon }_a{R}_L\downarrow -\epsilon R_L\uparrow$$

(8)

where $K\downarrow$ is the short-wave radiation. The $K\downarrow$ of the study region was observed by the Dongshan ground weather stations (seen Figure 1) and was approximately 680 Wm⁻². $R_L\downarrow$ and $R_L\uparrow$ are the downward and upward blackbody radiation in Wm⁻², respectively. $\alpha$ is the surface albedo, which was derived using the previously-mentioned conversion method proposed by Liang [38]. The resulting linear equation is:

$$\mathsf{\alpha }=0.356{\alpha }_1+0.130{\alpha }_3+0.373{\alpha }_4+0.085{\alpha }_5+0.072{\alpha }_7-0.0018$$

(9)

where $\mathrm{i}$ ($\mathrm{i}$ = 1–7) indicates the corresponding LANDSAT-5 TM band surface reflectance values.

$\epsilon$ is the surface emissivity. ${\epsilon }_a$ is the atmospheric emissivity, which can be calculated based on the following empirical equation proposed by [39]:

$${\epsilon }_a=1.24{(\frac{e_a}{T_a})}^{\frac{1}{7}}$$

(10)

where $e_a$ is the atmospheric water vapor pressure in hPa and $T_a$ is the atmospheric temperature in K. The short-wave radiation can be calculated based on sunshine duration according to [40]. The long-wave radiation can be calculated using Stefan–Boltzmann’s law:

$$R_L\downarrow =\sigma {T_a}^4$$

(11)

$$R_L\uparrow =\sigma {T_s}^4$$

(12)

where $\sigma$ is the Stefan–Boltzmann constant and $T_s$ is the surface temperature in K.

In this study, an advanced two-source energy model was used to calculate $H$ and $LE$. This method can decompose the heat fluxes within a mixed pixel into the sub-component heat fluxes of the vegetation and the non-vegetation.

$H$ is obtained using the bulk resistance method via the following equation:

$$H={\displaystyle \sum _{k=1}^m{F}_{veg}^k\times H_{veg}^k}+{\displaystyle \sum _{k=1}^n{F}_{non-veg}^k\times H_{non-veg}^k}$$

(13)

where $H_{veg}^k$ is the sensible heat flux for vegetated areas corresponding to vegetation endmember $k$ and $F_{veg}^k$ is the fraction of vegetation endmember $k$ within the pixel. $H_{non-veg}^k$ is the sensible heat flux for the non-vegetated areas corresponding to non-vegetation endmember $k$, and $F_{non-veg}^k$ is the fraction of the non-vegetated endmember $k$ within the pixel. $H_{veg}$ and $H_{non-veg}$ can be calculated using the following equations:

$$\begin{array}{l}{H}_{veg}=\rho C_p\frac{T_o-T_a}{R_{a\_veg}}\\ H_{non-veg}=\rho C_p\frac{T_o-T_a}{R_{a\_non-veg}+R_s}\end{array}$$

(14)

where $p$ is the air density, $C_p$ is the specific heat of air at constant pressure and $T_o$ is the surface aerodynamic temperature. Because $T_o$ is difficult to obtain, this study used remote sensing-derived surface temperature for $T_o$. $R_{a\_veg}$ and $R_{a_{non}-veg}$ are the aerodynamic resistances of the vegetated and non-vegetated areas, which were retrieved using the following equation:

$$R_a=\frac{\ln \left(\frac{Z_m-d}{Z_{0m}}\right)\times \ln \left(\frac{Z_h-d}{Z_{0h}}\right)}{k^{2}u}$$

(15)

where $Z_m$ is the height of the wind measurements, $Z_h$ is the height of the humidity measurements, $d$ is the zero plane displacement height, $Z_{0m}$ is the roughness length governing the momentum transfer and $Z_{0h}$ is the roughness length governing the transfer of heat and vapor. $k$ is von Karman’s constant, and $\mathrm{u}$ is the wind speed (m·s⁻¹) at a height of 2 m.

$R_s$ is the resistance to heat flow in the boundary layer immediately above the soil surface and can be computed by:

$$R_s=\frac{1}{a+b{u}_s}$$

(16)

where $a$ is the free convective velocity, $b$ is a coefficient that represents the typical soil surface roughness and $u_s$ is the wind speed over the soil surface at a height of 0.05–0.2 m. According to [15,41], $a$ and $b$ were set to 0.004 m/s and 0.012, respectively.

For the energy budget model used here, parameters involving the roughness lengths, $Z_{0m}$ and $Z_{0h}$, were important for the final results. Typical values of $Z_{0m}$ and $Z_{0h}$ for specific surface land types were alternatively used in this study. The surface types were obtained from image classification using the support vector machine method. Based on the associated existing literature [42,43,44], the parameters required in this study are given in Table 2.

Table 2. Parameters used for surface coverage types. LCT, land cover type.

LCTs	$Z_{0m}$	$Z_{0m}/Z_{0h}$	$d$
Water	0.3 × 10−4	0.34	0
Bare soil	0.001	50	0
Crop field	0.1	100	0.02
Lawn	0.001	50	0.133
Forest	0.5–1.0	1000	4
Developed, low	0.5	1000	2.5
Developed, high	1.5	1000	7.5

Table 2. Parameters used for surface coverage types. LCT, land cover type.

LCTs	$Z_{0m}$	$Z_{0m}/Z_{0h}$	$d$
Water	0.3 × 10−4	0.34	0
Bare soil	0.001	50	0
Crop field	0.1	100	0.02
Lawn	0.001	50	0.133
Forest	0.5–1.0	1000	4
Developed, low	0.5	1000	2.5
Developed, high	1.5	1000	7.5

The latent heat flux LE is calculated using:

$$LE={\displaystyle \sum _{k=1}^m{F}_{veg}^k\times L{E}_{veg}^k}+{\displaystyle \sum _{k=1}^n{F}_{non-veg}^k\times L{E}_{non-veg}^k}$$

(17)

where $L{E}_{veg}^k$ is the sensible heat flux for vegetated areas corresponding to vegetation endmember $k$ and $F_{veg}^k$ is the fraction of vegetation endmember $k$ within the pixel. $L{E}_{non-veg}^k$ is the sensible heat flux for the non-vegetated areas corresponding to vegetation endmember $k$, and $F_{non-veg}^k$ is the fraction of non-vegetated endmember $k$ within the pixel. $L{E}_{veg}$ and $L{E}_{non-veg}$ were computed as follows [41,45]:

$$\begin{array}{l}L{E}_{veg}=\frac{\rho C_p}{\gamma }\times \frac{e^o-e_a}{R_{a\_veg}+r_{s\_veg}}\\ L{E}_{non-veg}=\frac{\rho C_p}{\gamma }\times \frac{e^o-e_a}{R_{a\_non-veg}+r_{s\_non-veg}}\end{array}$$

(18)

where $e^o$ is the saturation vapor pressure in hPa and $\mathsf{\gamma }$ is the psychometric constant. $r_{s\_veg}$ and $r_{s{\_}_{non-veg}}$ are the stomata resistances for the vegetated areas and non-vegetated areas, respectively.

To obtain the stomata resistances $r_{s\_veg}$ and $r_{s{\_}_{non}-veg}$, this study used a common, but efficient method according to [14,17,18,46]. We assumed that the stomata resistance can be expressed by two environmental factors, air temperature and photo-synthetically active radiation, both of which can be retrieved from the remotely-sensed satellite data and ground meteorological data:

$$r_{s\_non-veg}=\frac{f_1(T_a)f_2(\text{PAR})}{r_{sMIN}}+\frac{1}{r_{cuticle}}$$

(19)

where PAR is the photosynthetic active radiation in Wm⁻², $r_{sMIN}$ is the minimum stomatal resistance in sm⁻¹ and $r_{cuticle}$ is the canopy resistance in sm⁻¹. $f_1$ and $f_2$ are retrieved by the equations developed in [47,48]. Following the work of [49] and [17], $r_{sMIN}$ was determined for each land cover type that was interpolated via the LCT classification of SVM.

In order to calculate the net heat storage on the urban surface, the storage heat flux (DG) was estimated by merging the ground heat flux and anthropogenic heat discharge based on the heat balance Equation (20). On the other hand, it should be noted that accurately discriminating the anthropogenic heat discharge from the ground heat flux is impractical because the quantitative differences between them are slight [8,18].

$$DG=G-A=Rn-H-LE$$

(20)

2.2.3. Retrieving the VFC, ISA% and Land Classification Map

MESMA was used to calculate the fractions of vegetation cover and impervious surfaces. Compared to other simple linear spectral mixture analyses that use the same types and number of endmembers to model all pixels, MESMA can decompose each pixel using different combinations of the represented endmembers. This approach was conducted in three stages. The pixel purity index method [50] was first used to find the most spectrally pure pixels (endmember) in the image used in this study. Once the purest pixels were chosen, we identified each selected endmember type based on their spectral characterization. A library containing 12 endmembers of crop fields, grasslands, forests, bare soils and impervious materials (dark and bright) was established in this study. The vegetation category included 6 distinct endmember spectra, while the bare soils and impervious surfaces included 2 and 4 spectra, respectively, as shown in Figure 3. Then, following Equations (21) and (22) [51], the collected endmembers were applied to the linear spectral unmixing analysis through an iterative process. The image was mixed by varying the material types of each class (i.e., vegetation and impervious surfaces) and the number of endmembers within each class. Moreover, the complexities of the combination models (i.e., two, three or four endmember models) were always varied until the fraction images achieved an accurate result. According to [52], three model constraints were imposed to choose the optimal resultant fractions for each pixel. Specifically, the models were considered acceptable if they met fraction and RMSE criteria, including that the fractions were physically reasonable (0–100% and added to 1) and that the RMSE was less than 2.5%.

$$R={\displaystyle \sum _{k=1}^n{f}_k{R}_k}$$

(21)

$${\displaystyle \sum _{k=1}^n{f}_k=1,}{f}_k\ge 0$$

(22)

where $k$ is the number of endmembers, $f_k$ is the fraction of endmember $k$ and $R_k$ is the reflectance of endmember $k$.

In this research, the VIPER Tools open-software package [53] was used to implement the MESMA algorithm. One hundred and four combination models were employed to characterize each pixel: 12 one-model combinations, 44 two-model combinations and 48 three-model combinations. Finally, an image with 12 bands that represented the fractions of each endmember was created, and VFC and ISA% were calculated as the sum of the fractions of the vegetation and impervious surfaces.

The fraction images derived from MESMA were combined with the classification result from the support vector machine (SVM) algorithm to classify the study region into seven classes: water, bare soils, crop fields, lawns, forests, low-medium density ISAs and high-density ISAs. Two types of impervious surfaces were distinguished considering the impervious surface density. According to a rule-based analysis, mixed pixels were classified as low-medium density ISAs when the percent ISA was less than 70, whereas those pixels greater than 70 were considered to be high-density ISAs.

Based on the above results, the two heat flux components (LE and H) were retrieved using a linear superposition of the representative sub-pixel land-cover heat fluxes. The two biophysical compositions (VFC and ISA%) were used here to compute the heat fluxes according to Equations (13) and (17).

Figure 3. Spectral reflectance of endmembers used for the MESMA.

Figure 3. Spectral reflectance of endmembers used for the MESMA.

3. Results and Discussion

3.1. Validation of the Land Cover Fractions and Heat Fluxes

To evaluate the accuracy of the predicted abundance images of vegetation cover and impervious surfaces using MESMA (see Figure 4a,b), following our previous studies [54,55], 100 sample plots were randomly selected, and the corresponding reference values were obtained from a Google Earth™ image from 12 May 2010.

Figure 4. Maps of land coverage, LCTs classification and LCT at 1040 LT, derived from 30-m LANDSAT-5 TM acquired on 24, May 2010. (a) Vegetation fraction coverage; (b) impervious surface coverage; (c) LCT classification; (d) land surface temperature.

Figure 4. Maps of land coverage, LCTs classification and LCT at 1040 LT, derived from 30-m LANDSAT-5 TM acquired on 24, May 2010. (a) Vegetation fraction coverage; (b) impervious surface coverage; (c) LCT classification; (d) land surface temperature.

As illustrated in Figure 5, the regression analysis showed a relationship with an R² of 0.85 and a slope of 0.96 regarding the VFC. Concerning the VFC, the overall RMSE was 7.79, and the overall bias was 6.40. RMSEs of 7.16 and 9.55 were achieved for the high-density impervious surface areas and the low-medium density impervious surface areas, respectively. Biases of 5.87 and 8.21 were achieved for high-density impervious surface areas and the low-medium density impervious surface areas, respectively. For the ISA%, the R² was 0.83, and the slope was 0.86. The results exhibited an overall RMSE of 8.67, and an overall bias of 6.8. RMSE of 7.41 and 12.34 were achieved for the low-medium density impervious surface areas and high-density impervious surface areas, respectively. Biases of 5.59 and 10.80 were achieved for the low-medium density impervious surface areas and the high-density impervious surface areas, respectively. Compared to the results of our previous studies that were based on simple spectral mixture analyses [54,55], the results achieved in this study were acceptable.

Figure 5. Accuracy assessment results of (a) estimated vegetation fraction coverage (VFC) and (b) estimated impervious surface coverage (ISA%).

Figure 5. Accuracy assessment results of (a) estimated vegetation fraction coverage (VFC) and (b) estimated impervious surface coverage (ISA%).

The LCT classification obtained using MESMA and SVM is shown in Figure 4c. The LST retrieved using the LANDSAT-5 TM image of 24 May 2010 is shown in Figure 4d. It was found that although the identification of some dominant impervious surfaces of complex urban landscapes was relatively less successful, distinctions among most of the other classes, including vegetation cover, bare soils and water bodies, were observed. Specifically, the SVM method yielded an overall accuracy of 90.1% and a kappa value of 0.88. The producer’s accuracies for crops, grasslands and forests were 87%, 89% and 94%, respectively. Overall, the combination of SVM and MESMA was able to achieve reasonably sound LCT results in the present study.

Figure 6 shows the surface heat fluxes from the LANDSAT-5 TM image of 24 May 2010. Because there was no in situ eddy flux measurement over Suzhou in 2010, unfortunately, two alternative methods were adopted to evaluate the calculated results. The first method compared the estimated daily ET with pan evaporation measurements taken at meteorological stations following the method of [56] and [57]. The instantaneous LEs at two stations (Suzhou and Dongshan) were extended to daily ETs using the instantaneous evaporative fraction and 24-h net radiation. The calculated daily ETs were 3.89 mm·day⁻¹ at the Suzhou station and 4.13 mm·day⁻¹ at the Dongshan station, and the daily ETs measured using the evaporation pan at the Suzhou and Dongshan stations were both 4.0 mm·day⁻¹. It should be kept in mind that daily ETs collected from pan evaporation measurements are different from the actual ETs and are largely affected by factors that include local landscape position and pan maintenance. Given this, the estimated daily ET was also assessed using another relative qualitative criterion. The relative errors between the measured and retrieved values were 2.75% and −3.25%, respectively, which indicates that our method can produce an acceptance flux result. However, in the future, the predicated daily ETs need to be validated using more reliable measurements, such as the ET measured from the eddy-covariance technique. Meanwhile, the resultant heat fluxes were compared to previously published findings at other sites in terms of the ratios of turbulent heat fluxes. Table 3 shows that the mean H/Rn ratio was approximately 0.20 in the urban areas, and the mean LE/Rn ratios were 0.42 and 0.1 in the high density ISAs and low-medium density ISAs, respectively. The tendencies for both H/Rn and LE/Rn are basically similar to those in [8,13,14], which implies that the results are reasonable. It should be mentioned that because LE is the energy used to move water vapor from land surfaces to the atmosphere through vegetation evapotranspiration or soil evaporation, it is understandable that the mean ratio of LE/Rn in Suzhou is somewhat higher in other studies, such as [58].

In addition, H was recorded to be approximately 20% of the Rn in the urban regions, but below 10% of the Rn over the rural areas. The mean ratio of H/Rn for the urban areas was lower in comparison with previous studies, probably due to differences in structural and material properties of buildings. The mean ratios of DG/Rn were 0.47 and 0.68 for low-medium density ISAs and high-density ISAs, respectively. These observations were slightly inconsistent with those of other studies, probably because anthropogenic heating is significantly affected by factors such as landscape, albedo and building areas [59].

Figure 6. The distribution maps of heat fluxes (in Wm⁻²) at 1040 LT, estimated with 30-m LANDSAT-5 TM acquired on 24 May 2010. (a) Net radiation, Rn; (b) sensible heat flux, H; (c) latent heat flux, LE; and (d) storage heat flux, DG.

Figure 6. The distribution maps of heat fluxes (in Wm⁻²) at 1040 LT, estimated with 30-m LANDSAT-5 TM acquired on 24 May 2010. (a) Net radiation, Rn; (b) sensible heat flux, H; (c) latent heat flux, LE; and (d) storage heat flux, DG.

Table 3. Mean ratios of heat fluxes to net radiation over different LCTs.

Mean Value	H/Rn (%)	LE/Rn (%)	DG/Rn (%)
Bare soil	3.81	32.5	63.81
Crop	5.56	65.6	33.22
Lawn	2.35	53.1	44.58
Forest	3.60	85.85	15.27
Low-medium density ISA	18.78	41.38	46.71
High density ISA	22.38	10.0	68.07

Table 3. Mean ratios of heat fluxes to net radiation over different LCTs.

Mean Value	H/Rn (%)	LE/Rn (%)	DG/Rn (%)
Bare soil	3.81	32.5	63.81
Crop	5.56	65.6	33.22
Lawn	2.35	53.1	44.58
Forest	3.60	85.85	15.27
Low-medium density ISA	18.78	41.38	46.71
High density ISA	22.38	10.0	68.07

3.2. Spatial Pattern of Heat Energy

The mean values and standard deviations (SDs) of the LST and heat fluxes (Rn, LE, H and DG) for each surface type were computed and are shown in Table 4. An examination of these values shows that a distinct thermal pattern was occurring among the different LCTs.

The spatial pattern of LST over Suzhou city is shown in Figure 5d, which shows that the urban areas captured the highest LST. The LST tended to decrease from downtown to the surrounding suburban areas. The mean LST in the urban areas was approximately 4 K higher than that in the suburban areas, which indicates that an obvious UHI effect occurred in Suzhou at 10:40 a.m. on 24 May 2010. The vegetation cover had a lower mean value and smaller variations in LST due to its thermal capacity for reducing heat storage. Moreover, the mean LST of the high density impervious surface areas was larger than that of the low-medium density impervious surface areas. The LSTs of the former had a mean value of 303.3 K and an SD of 2.0 K, whereas the mean value and SD of the latter were 302.6 K and 1.8 K, respectively. This difference could be attributed to the difference in LCTs compositions. Another cause that cannot be ignored is the differences in dense populations over different developed urban surfaces [60,61].

Table 4. Descriptive statistics for heat fluxes over different LCTs.

	Rn		H		LE		DG		LST
	Mean	SD	Mean	SD	Mean	SD	Mean	SD	Mean	SD
Water	593.5	18.5	−8.9	7.6	477.3	29.2	118.4	16.2	295.6	2.1
Bare soil	490.2	35.6	17.8	8.6	157.1	31.4	315.5	68.3	302.5	2.1
Crop	528.8	28.1	21.1	23.8	352.2	72.3	169.4	111.6	299.8	2.1
Lawn	474.3	17.3	10.9	5.8	282.8	19.8	190.6	361	300.8	1.5
Forest	534.2	19.2	13.2	19.5	462.6	46.9	59.4	73.6	298.8	1.2
Low-medium density ISA	491.5	24.2	90.4	36.7	207.9	85.7	214.1	95.8	302.6	1.8
High density ISA	485.7	46.4	104.9	40.2	47.8	72.7	334.0	106.5	303.3	2.0

Table 4. Descriptive statistics for heat fluxes over different LCTs.

	Rn		H		LE		DG		LST
	Mean	SD	Mean	SD	Mean	SD	Mean	SD	Mean	SD
Water	593.5	18.5	−8.9	7.6	477.3	29.2	118.4	16.2	295.6	2.1
Bare soil	490.2	35.6	17.8	8.6	157.1	31.4	315.5	68.3	302.5	2.1
Crop	528.8	28.1	21.1	23.8	352.2	72.3	169.4	111.6	299.8	2.1
Lawn	474.3	17.3	10.9	5.8	282.8	19.8	190.6	361	300.8	1.5
Forest	534.2	19.2	13.2	19.5	462.6	46.9	59.4	73.6	298.8	1.2
Low-medium density ISA	491.5	24.2	90.4	36.7	207.9	85.7	214.1	95.8	302.6	1.8
High density ISA	485.7	46.4	104.9	40.2	47.8	72.7	334.0	106.5	303.3	2.0

The Rn was found to be higher for the vegetation covers (crops, grass and forests) and was relatively lower for the impervious surface areas (Figure 6a). Furthermore, a high Rn was found in the forest cover because trees among all vegetation types can absorb the most solar radiation. Low Rns were observed in the bare soil and impervious surfaces, whose land cover types are generally accompanied by high albedos. In addition, water pools were found to have high Rns, with a mean value of 593.47 Wm⁻² due to their low LSTs and surface albedos.

An examination of the mean value of H found that impervious surfaces had higher Hs than those of vegetation cover. Meanwhile, the SD of H for impervious surfaces was higher than that of vegetation cover. Furthermore, the high density ISAs had obviously higher H relative to those of the low-medium density ISAs, which indicates that more energy was redistributed between the surface and the atmospheric boundary through convection in highly developed land. The differences in H among the different density impervious areas may provide a possible way to mitigate the UHI effect by changing the land compositions of urban landscapes.

Concerning LE, its mean values for the crop fields, grasslands, forests and water were high. LE was found to be lower for bare land and impervious areas, which both had low water availability. It is not surprising that the city of Suzhou has a smaller mean LE in the downtown urban area than in the surrounding rural areas, because the LE is mainly controlled by vegetation cover rather than impervious surfaces. Moreover, although the mean LE was lower over the impervious surface areas due to the decreased vegetation cover percentage, it still exhibited a higher mean value in urban green areas.

Because the fraction of vegetation cover distributed in Suzhou will perform evapotranspiration, the VFC retrieved from MESMA had to be used to estimate the LE for the pixels that contained mostly, but not exclusively impervious surfaces. Table 4 shows that the mean value of LE in low-medium density impervious surface areas was 207.9 Wm⁻². Moreover, unlike most remote sensing studies that assigned an unreasonable value of zero to the mixed pixels, the high density impervious surface areas in the present study had non-zero values of LE (mean value of 47.8 Wm⁻²).

Regarding the storage heat flux, the increased intensity of which was considered by [62] to be a significant feature for urban energy processes, its mean values were 334.0 Wm⁻² and 214.1 Wm⁻² in the high-density ISAs and low-medium density ISAs, respectively. The mean values of DG for crops, grasslands and forests were 169.4 Wm⁻², 190.6 Wm⁻² and 59.4 Wm⁻², respectively. Moreover, it has been observed that impervious surfaces can store more heat than the surrounding vegetation cover. This difference can be explained by the fact that impervious surfaces and vegetation cover are two different types of land cover with distinct thermal properties, such as thermal conductivity and heat capacity. On the other hand, [8,10] reported that the surface enlargement in built up areas caused by the three-dimensional fabric can directly add urban areas and volume to the shallow depth layers, consequently leading to the storage of more energy. In addition, the regions of beach land mostly were found to have a higher DG value. This is related to the fact that surface DG is closely associated with soil moisture. Furthermore, despite the higher DG in some industrial areas, the DG was relatively low at some inland factories. It is likely that low DG in industrial areas is the result of energy consumption.

3.3. The Relationships between Heat Energy and VFC and ISA%

Different categories of VFC and ISA% were computed to analyze the relationships among VFC, ISA% and heat energy (LST and heat fluxes). The mean values and SDs of LST, LE and H for each VFC and ISA% category are shown in Table 5 and Table 6. In addition, integrating the urban heat energy with VFC and ISA% with appropriate intervals may yield interesting findings. Therefore, on the basis of prior studies that examined the urban thermal environment [54,63,64], a zonal analysis was performed to examine the mean value of LST and heat flux on an increment of vegetation coverage from 0.1–1 and with an increment of ISA% from 10–100.

Table 5. Mean LE and LST of vegetation coverage categories.

	LE		LST
VFC	Mean	SD	Mean	SD
<10%	66.5	122.8	302.6	2.7
10%–30%	174.2	122.1	302.6	2.1
30%–50%	253.3	109.7	302.5	2.3
50%–70%	325.3	98.6	301.3	2.3
>70%	361.3	81.7	299.9	2.1

Table 5. Mean LE and LST of vegetation coverage categories.

	LE		LST
VFC	Mean	SD	Mean	SD
<10%	66.5	122.8	302.6	2.7
10%–30%	174.2	122.1	302.6	2.1
30%–50%	253.3	109.7	302.5	2.3
50%–70%	325.3	98.6	301.3	2.3
>70%	361.3	81.7	299.9	2.1

Table 6. Mean H and LST of percentage of ISA categories.

	H		LST
ISA	Mean	SD	Mean	SD
<10%	15.8	23.1	299.4	1.9
10%–30%	44.2	40.1	300.9	2.3
30%–50%	61.1	46.0	301.8	2.2
50%–70%	68.0	49.6	302.8	2.1
>70%	88.1	48.9	303.0	2.2

Table 6. Mean H and LST of percentage of ISA categories.

	H		LST
ISA	Mean	SD	Mean	SD
<10%	15.8	23.1	299.4	1.9
10%–30%	44.2	40.1	300.9	2.3
30%–50%	61.1	46.0	301.8	2.2
50%–70%	68.0	49.6	302.8	2.1
>70%	88.1	48.9	303.0	2.2

The correlations between LST and VFC and ISA% are illustrated in Figure 7. The results show that the LST value was negatively correlated with VFC, but was positively correlated with ISA%. A statistically-significant correlation (p < 0.05; adjusted $R^2$ = 0.92) was found between the mean LST and VFC, and a statistically-significant correlation (p < 0.05; adjusted $R^2$ = 0.84) was also found between the mean LST and ISA%. These observations suggest that both the vegetation cover and impervious surface substantially contribute to urban surface energy exchange.

Figure 7. Scatterplots of (a) mean LST versus vegetation coverage and (b) mean LST versus percentage of ISA.

Figure 7. Scatterplots of (a) mean LST versus vegetation coverage and (b) mean LST versus percentage of ISA.

Table 5 reveals that the high-density vegetation cover had a higher mean LE than those of the other categories. Figure 7 shows a scatterplot of the mean LE versus the corresponding mean VFC that illustrates the statistically-significant correlation (p < 0.05; adjusted $R^2$ = 0.88) between them. This strong positive correlation indicates that the variations in latent heat flux can be well delineated by urban vegetation coverage. On the other hand, we found that the high-density impervious urban areas had higher mean values of H relative to those of the lower density zones. Figure 8 shows a scatterplot of the mean H and corresponding ISA%, and the result indicates that there was a statistically-significant correlation (p < 0.05; adjusted $R^2$ = 0.90) between them. In the study of [10], the relationships between the heat fluxes and VFC and ISA% were characterized in terms of the linearization of the Bowen ratio. The present study shows that urban green areas can effectively regulate urban surface temperatures, not only by absorbing solar radiation, but also through the cooling effects of latent heat [65,66]. However, the correlations among them may be subject to many factors [67], and more attention should therefore be given to combined field surveys and computer simulation modeling in future studies.

In addition, it can be seen that the VFC-LE relationship tended to be more significant for the sparse vegetation cover than for the dense vegetation cover. Figure 8 reveals that vegetation cover was more substantially correlated with LE in the low-density vegetation regions, which indicates that lower vegetation coverage can lead to a more efficient impact on the surface LE. This may be ascribed to advection effects and the serious shading caused by various water pools and large buildings during the daytime. Similar phenomena were also found by the studies of [68,69], which reported that green spaces distributed in intensively built-up areas were prone to generate high latent heat due to microscale advection. In addition, for the impervious surface areas around the water bodies, which occupied a considerable fraction of the study region, it was also found that the surface DG was reduced due to the increased H by the advection of cooler air from the water bodies. Therefore, municipal governments should focus on increasing green areas and small water bodies in city centers, especially for intensively built-up lands, such as central business districts (CBDs).

Figure 8. Scatterplots of (a) mean LE versus vegetation coverage and (b) mean H versus percentage of ISA.

Figure 8. Scatterplots of (a) mean LE versus vegetation coverage and (b) mean H versus percentage of ISA.

3.4. Significance and Future Work

The LE of impervious surfaces retrieved in most remote sensing studies was set to zero based on the assumption that the impervious materials were dry when the satellite images were acquired. Moreover, the approach based on traditional hard classification schemes would cause an underestimation of LE, especially for pixels that contain predominately, but not exclusively impervious surfaces. From this perspective, the present study attempted to obtain LE and H using a method that considers per-pixel heat fluxes to be a linear combination of the representative sub-pixel heat fluxes.

On the other hand, although fewer studies have used the MESMA approach to calculate urban heat fluxes, simple spectral mixture analyses have proven to be efficient for describing urban heat flux terms [70,71]. Qualitatively, our results complement the findings of a small number of previous studies that revealed the potential of spectral mixture analysis for depicting urban surface energy fluxes. However, this approach is sensitive to the collection of representative spectral data. Endmember selection generally determines how accurately mixture analysis models can represent spectra [72,73]. In this paper, the 12 most optimal spectra were obtained from the utilized image in consideration of the operational efficiency. Consistent with the studies of [74,75], we believe that the performance of MESMA and further modeling analysis will be enhanced if more representative spectral data, especially of urban impervious surfaces and green areas, are used, and this will be studied in future studies.

The aerodynamic roughness parameters of urban LCTs are essential for describing the energy processes that occur within urban air boundary layers and urban climates. Unfortunately, these parameters (the roughness length and zero-plane displacement) were solely obtained from a look-up table in our work and thus may have affected the accuracies of the predicated LE and H. The work in [17,18] comprehensively assessed the importance of aerodynamic roughness for estimating urban heat flux. Accordingly, we believe that both H and LE can be more accurately retrieved using the bulk resistance method when the roughness length is acquired by field observations or with higher spatial resolution remote sensing data, which deserves further consideration.

The LANDSAT-5 TM image, which has spatial resolutions of 30 m and 120 m for the optical and thermal bands, respectively, was used to retrieve the LCT classification, LST and heat fluxes. The moderate spatial resolution (especially 120 m for LST) may not be fine enough for quantitatively examining urban surface energy fluxes. However, the overall accuracy can be improved if complex urban landscapes are classified using higher spatial resolution images or the LST is derived from finer spatial resolution thermal images. Some new types of sensors (i.e., GaoFen and Sentinel), from which more abundant and new types of high spatial resolution remote sensing data can be acquired, will help solve this limitation in the future.

4. Conclusions

This study estimated urban heat fluxes using an advanced two-source method. In this method, the heat fluxes of sub-pixel vegetation covered areas were taken into account, and the heat fluxes (LE and H) of the non-homogeneous pixels were regarded as a combination of the sub-pixel level heat fluxes and their corresponding land-cover area proportions. The accuracies of the resulting heat fluxes were assessed using evaporation pan data from meteorological stations and the mean ratios of the turbulent heat fluxes and the ground storage heat (DG) to the net radiation (Rn). The results show that the estimated heat fluxes were plausible and can be used for further urban thermal environment analyses.

Our study also has enhanced the understanding of the spatial patterns of urban heat energy in response to different LCTs and the two most common urban biophysical compositions. It should be concluded that the environmental variables, VFC, ISA% and LCT, are promising for quantifying patterns of urban heat flux.