Simulating the Impact of Urban Surface Evapotranspiration on the Urban Heat Island E ﬀ ect Using the Modiﬁed RS-PM Model: A Case Study of Xuzhou, China

: As an important energy absorption process in the Earth’s surface energy balance, evapotranspiration (ET) from vegetation and bare soil plays an important role in regulating the environmental temperatures. However, little research has been done to explore the cooling e ﬀ ect of ET on the urban heat island (UHI) due to the lack of appropriate remote-sensing-based estimation models for complex urban surface. Here, we apply the modiﬁed remote sensing Penman–Monteith (RS-PM) model (also known as the urban RS-PM model), which has provided a new regional ET estimation method with the better accuracy for the urban complex underlying surface. Focusing on the city of Xuzhou in China, ET and land surface temperature (LST) were inversed by using 10 Landsat 8 images during 2014–2018. The impact of ET on LST was then analyzed and quantiﬁed through statistical and spatial analyses. The results indicate that: (1) The alleviating e ﬀ ect of ET on the UHI was stronger during the warmest months of the year (May–October) but not during the colder months (November–March); (2) ET had the most signiﬁcant alleviating e ﬀ ect on the UHI e ﬀ ect in those regions with the highest ET intensities; and (3) in regions with high ET intensities and their surrounding areas (within a radius of 150 m), variation in ET was a key factor for UHI regulation; a 10 W · m − 2 increase in ET equated to 0.56 K decrease in LST. These ﬁndings provide a new perspective for the improvement of urban thermal comfort, which can be applied to urban management, planning, and natural design.


Introduction
The urban heat island (UHI) effect describes the phenomenon through which the urban temperatures are significantly higher than in the surrounding suburbs [1]. Urbanization is typically accompanied by the transformation of vegetation, soil, and other surface cover types into cement, asphalt, concrete, and other impervious surfaces [2]. Widespread changes in surface and the emission of anthropogenic heat in urban areas not only change regional surface radiation and heat conduction characteristics but also have a significant impact on the regional energy balance, water cycle, and atmospheric conditions. In addition, urban population density [3][4][5], urban energy consumption [6], urban climatic and meteorological conditions [5,7], and urban atmospheric pollution [8] have been widely shown to be the main factors influencing the formation and enhancement of UHI effect. energy balance; and (3) the component surface characteristics parameters are added in the inversion of component net radiation, which has simulated the differences of the net radiations for various surfaces. Compared with the MMP model, the urban RS-PM model has greatly simplified the inversion process and provides an acceptable level of accuracy given that modeled ET values are less sensitive to variations in the component temperatures, which can be directly estimated using empirical formula. Furthermore, the urban RS-PM model can be applied to more types of remote sensing data, such as the Landsat series and the HJ-1 A/B satellites (manufactured by China Aerospace Science and Technology Corporation, and launched on September 6, 2008) [35]. Therefore, the urban RS-PM model offers an appropriate tool to explore the mechanisms through which ET mitigates the UHI at a regional scale.

Study Area
Xuzhou is located between 117°01'-117°25'E and 34°06'-34°22N ( Figure 1) and is the secondlargest city in the northwestern province of Jiangsu, East China. According to 2018 statistics, Xuzhou covers a total area of 11,258 km 2 , including a 3037-km 2 urban area, with a population of 87.63 million of which 55.88 million live in urban areas [40]. Xuzhou has a warm temperate monsoon climate with annual precipitation between 800 and 930 mm, and a rainy season that accounts for approximately 56% of the annual total precipitation. The annual average percentage of sunshine is between 52% and 57%, and the annual average temperature is approximately 15.5 °C (2014-2016 data). Vegetation cover is approximately 30.3%, and the urban area coverage has reached 43.3%.
As a typical medium-sized city in China, Xuzhou has undergone large-scale urbanization but has relatively good environmental (climatic and ecological) conditions due to its rich vegetation coverage that is conducive to surface ET.  As a typical medium-sized city in China, Xuzhou has undergone large-scale urbanization but has relatively good environmental (climatic and ecological) conditions due to its rich vegetation coverage that is conducive to surface ET.

Satellite Data
Ten sets of Landsat 8 Operational Land Imager (OLI, 30 m resolution) and Thermal Infrared Sensor (TIRS, 100 m resolution) images (all of the satellite datasets used in our study are listed in Table 1) for the period 2014-2018 were selected for inversing urban surface evapotranspiration and heat island effect in the study area. Six image sets were acquired during (or close to) summer (May-October, classified as warm season) and four during (or close to) winter (November-March, classified as cold season). For the Landsat 8 data, OLI bands 1-7 were used to extract the pixel component fractions and the intermediate parameters in the ET inversion, and TIRS Band 10 was used to retrieve land surface temperatures. Based on the acquisition date of the Landsat 8 data, four GF-1 and GF-2 high-resolution remote sensing images were also selected to validate the accuracy of the linear spectral mixture analysis. Based on the urban RS-PM model, it is necessary to obtain meteorological observation data including air temperature (T air ), air relative humidity (RH), atmospheric pressure (P A ), wind speed (u z ), and water vapor pressure (e) for inversing urban surface ET. All meteorological data (Table 2) corresponding to the acquisition time of the Landsat 8 data were obtained from the meteorological observatory at the China University of Mining and Technology (CUMT) in Xuzhou, and the time resolution of the meteorological data is 30 min.

Flux Observations
To validate the accuracy of the ET values estimated using the urban RS-PM model, observed latent heat flux data were also obtained from the open-path eddy covariance (EC) system set up on a 30 m-high flux tower located at the CUMT Collaborative Observation Test Site ( Figure 2). According to the aerodynamic environment around the flux tower, the height above the ground (Z) of the EC was 15 m, the average building height around the flux tower (h b ) was 4.5 m, the zero displacement (d) was 0.68 m, the momentum roughness length (Z om ) was 3 m, and the effective measurement height (Z − d) was 14.32 m. To validate the accuracy of the ET values estimated using the urban RS-PM model, observed latent heat flux data were also obtained from the open-path eddy covariance (EC) system set up on a 30 m-high flux tower located at the CUMT Collaborative Observation Test Site ( Figure 2). According to the aerodynamic environment around the flux tower, the height above the ground ( ) of the EC was 15 m, the average building height around the flux tower (ℎ ) was 4.5 m, the zero displacement ( ) was 0.68 m, the momentum roughness length ( ) was 3 m, and the effective measurement height ( − ) was 14.32 m.

ET Estimation Using the Urban RS-PM Model
The urban RS-PM model [35] was developed from the surface energy balance dual-source parallel model [33] and the RS-PM model [26]. Mixed pixels of Landsat 8 images of urban areas are composed of impervious surface, vegetation, and soil components, while water generally exists in isolated patches and can be masked separately [41]. As an important component of the urban mixed pixels, impervious surface is an integral endmember in the linear spectral analysis for extracting vegetation and bare soil fractions. However, it cannot retain water in dry weather, which makes no contribution to urban evapotranspiration. Therefore, the ET from an urban mixed pixel can be estimated as the vegetation component transpiration and the soil component water evaporation, as follows:

ET Estimation Using the Urban RS-PM Model
The urban RS-PM model [35] was developed from the surface energy balance dual-source parallel model [33] and the RS-PM model [26]. Mixed pixels of Landsat 8 images of urban areas are composed of impervious surface, vegetation, and soil components, while water generally exists in isolated patches and can be masked separately [41]. As an important component of the urban mixed pixels, impervious surface is an integral endmember in the linear spectral analysis for extracting vegetation and bare soil fractions. However, it cannot retain water in dry weather, which makes no contribution to urban evapotranspiration. Therefore, the ET from an urban mixed pixel can be estimated as the vegetation component transpiration and the soil component water evaporation, as follows: Remote Sens. 2020, 12, 578 6 of 22 where ET v and ET s are the vegetation component transpiration and soil component evaporation in an urban mixed pixel, respectively; f v and f s are the vegetation and soil fractions in the mixed pixel, respectively; ET * v is transpiration in a pure vegetation pixel; and ET * s is evaporation in a pure soil pixel. A detailed explanation for Equations (1) to (3) is shown in Figure 3. and are the vegetation and soil fractions in the mixed pixel, respectively; * is transpiration in a pure vegetation pixel; and * is evaporation in a pure soil pixel. A detailed explanation for Equations (1) to (3) is shown in Figure 3. Considering the heterogeneity of pixel components, the urban RS-PM model used a re-derived Penman-Monteith algorithm that incorporates the surface energy balance dual-source parallel model. As such, the vegetation component transpiration and soil water evaporation ( and ) in an urban mixed pixel can be estimated as follows [35]: where (J·kg -1 ) is the latent heat of vaporization; ∆ is the slope of the relationship between saturated water vapor pressure ( , Pa) and temperature; , * and , * (W·m -2 ) are the net radiation in a pure vegetation pixel and a pure soil pixel, respectively; * (W·m -2 ) is the soil heat flux in a pure soil pixel; (kg·m -3 ), (J·kg -1 ·K -1 ), and (Pa) are the air density, air specific heat capacity, and actual water vapor pressure, respectively; , and , (s·m -1 ) are the aerodynamic resistance to heat transfer of vegetation and soil, respectively; is the psychrometric constant; , (Pa·K -1 ) is the vegetation canopy surface resistance; (Pa·K -1 ) is the total aerodynamic resistance to vapor transport at the soil surface (including the surface resistance and aerodynamic resistance for vapor transport); and is the air relative humidity. In Equations (4) and (5), the parameters , ∆, , , , , and can be directly obtained or calculated from meteorological observations. The other parameters ( , , Considering the heterogeneity of pixel components, the urban RS-PM model used a re-derived Penman-Monteith algorithm that incorporates the surface energy balance dual-source parallel model. As such, the vegetation component transpiration and soil water evaporation (ET v and ET s ) in an urban mixed pixel can be estimated as follows [35]: where λ (J·kg −1 ) is the latent heat of vaporization; ∆ is the slope of the relationship between saturated water vapor pressure (e s , Pa) and temperature; R * n,v and R * n,s (W·m −2 ) are the net radiation in a pure vegetation pixel and a pure soil pixel, respectively; G * s (W·m −2 ) is the soil heat flux in a pure soil pixel; ρ (kg·m −3 ), C p (J·kg −1 ·K −1 ), and e a (Pa) are the air density, air specific heat capacity, and actual water vapor pressure, respectively; r ah,v and r ah,s (s·m −1 ) are the aerodynamic resistance to heat transfer of vegetation and soil, respectively; γ is the psychrometric constant; r s,v (Pa·K −1 ) is the vegetation canopy surface resistance; r tot (Pa·K −1 ) is the total aerodynamic resistance to vapor transport at the soil surface (including the surface resistance and aerodynamic resistance for vapor transport); and RH is the air relative humidity.
In Equations (4) and (5), the parameters λ, ∆, ρ, C p , e s , e a , and RH can be directly obtained or calculated from meteorological observations. The other parameters ( f v , f s , R * n,v , R * n,s , r ah,v , r ah,s , r s,v , r s,s , and G s ) were calculated using the surface information extracted from remote sensing imagery.

Linear Spectral Analysis
A normalized linear spectral mixture analysis (NSMA) was applied (see Wu et al. [42]) to extract the proportions of vegetation ( f v ), soil ( f s ), high-albedo impervious surface ( f imp_h ), and low-albedo impervious surface ( f imp_l ) for each mixed pixel, which can be expressed as follows: where R b is the normalized reflectance of band b for a pixel; f i is the fraction of each endmember (where f 1 , f 2 , f 3 , and f 4 refer to f v , f s , f imp_h , and f imp_l , respectively); R i,b is the endmember normalized reflectance of band b for the pixel; N is the total number of all endmembers; and e b is the residual. It is necessary to use GF-1/2 high-resolution images to verify the accuracy of linear spectral analysis. One hundred sample points were selected randomly in the Landsat 8 image, and each sample point contained 3 × 3 pixels (90 × 90 m). Then, the GF-1/2 images were manually interpreted to extract the true component fractions of each sample point to verify the accuracy of the inversed f v and f s .

Component Net Radiation Inversion
The net radiation flux was determined as the incoming net shortwave radiation and net longwave radiation [7]. For a pure pixel of vegetation or soil, net radiation flux (R * n,v and R * n,s ) were calculated as follows [35,36,43]: where a v and a s are the typical surface albedos of pure vegetation and pure soil pixels, respectively, which can be extracted from Landsat 8 images [44]. Here, taking the average value of a v = 0.18 and a s = 0.28 for the study area, S d (W·m −2 ) is the incoming solar radiation; ε air is the atmospheric emissivity; σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W·m −2 ·K −4 ); T air (K) is the air temperature, which was estimated using Qin et al.'s [45] algorithm; ε v and ε s are the typical surface emissivity of vegetation and soil, respectively; and T v and T s (K) are the component surface temperature of vegetation and soil, respectively.
• Calculation of ε air In cloudless weather, the atmosphere emissivity ε air can be calculated by using the atmospheric water vapor pressure e and the atmospheric average temperature T air [50]: ε air = 1.24(e/T air ) 1/7 (19) • Estimation of T v and T s As there is only one TIR band for Landsat 8 data, an empirical formula of pixels' component radiation ratio [35,45,51] was used here to estimate component temperatures T v and T s : The method for inversing the land surface temperature T sur is shown in Section 2.4, and NDVI v and NDVI s are the NDVI threshold values of pure soil and pure vegetation pixels, respectively.

Component Aerodynamic Resistance Inversion
The algorithm for estimating aerodynamic resistance to heat transfer was developed by Kustas and Norman [52]. For vegetated and soil surfaces, aerodynamic resistance (r ah,v and r ah,s ) can be expressed as follows: where r ah,i is the aerodynamic resistance, and i = 1 and i = 2 refer to vegetation and soil, respectively; k is the von Karman constant (0.4); Z is the reference height; u z is the wind speed observed at height Z; d o,i is the zero-plane displacement; Z om,i and Z oh,i are the roughnesses for momentum and sensible heat, respectively; and ψ m and ψ h are the stability correction functions for momentum and heat, respectively, which were calculated from the Monin-Obukhov length L m [53].
• Estimation of ψ m,i For unstable conditions (∀Z/L M < 0): For stable conditions (∀Z/L M = 0): Brutsaert et al. [54] have proposed the algorithm for estimating Z om,v and d o,v by using vegetation height h v , which was applicable to most vegetation types [55]. Z oh,v can be estimated by Kustas et al.'s [55] algorithm.
where S kb is the empirical coefficient with S kb = 0.13; u z is the wind speed (m·s −1 ); T v is vegetation component temperature (K); and T air is atmospheric average temperature.
where Re * is the roughness Reynolds number; u * is the friction velocity; and v is the kinematic molecular viscosity with v = 1.48 × 10 −7 m 2 ·s −1 .

Component Surface Resistance Inversion
• Estimation of r s,v By calculating canopy conductance, Mu et al. [26] developed a method for estimating vegetation canopy surface resistance (r s,v ), as follows: where the leaf area index (LAI) was used to convert stomatal conductance (C s ) to canopy conductance (C c ). The Level-4 MODIS global leaf area index product (MOD15A2H) was used for LAI data, which were acquired from Level 1 and Atmosphere Archive and Distribution System Distributed Active Archive Center (LAADS DAAC) [59]; m(T min ) is a constraint that limits C s via minimum air temperatures (T min ); m(VPD) is the constraint for reducing C s relative to the saturated air vapor pressure (VPD); and C L is the mean potential stomatal conductance per unit leaf area which can be taken as 0.0013 according to Leuning et al.'s [60] research.
where VPD open , VPD close , T min_open , and T min_close are the threshold values under the conditions of nearly complete suppression with full stomatal closure ("close" condition) and no suppression for transpiration ("open" condition), respectively. The threshold values for various biomes can be obtained from the Biome Properties Look-Up Table (BPLUT) [32].
• Estimation of r tot For soil surfaces, Mu et al. [26] proposed that the total aerodynamic resistance for vapor transport (r tot ) can be used to estimate soil transport. They assumed a global constant for r tot based on ground observations [61] and corrected this using standardized conditions; thus: r totc = 107.0 m·s −1 (40) r tot = r totc × r corr (41) where r corr is the correction term; r totc is the constant based on ground observations; and T ( • C) and P A (kPa) are the atmospheric temperature and pressure, respectively.

Soil Heat Flux Inversion
Based on the observation data of the soil heat flux and the soil component net radiation, Friedl [62] has reported that there is a functional relationship between soil heat flux G * s , soil component net radiation R * n,s , and solar zenith angle θ, which can be expressed as follows:

Land Surface Temperature Inversion Using the IMW Algorithm
As there are two TIR bands (Band 10 and Band 11) in Landsat 8 data, the split-window algorithm is suitable for inversing land surface temperature (LST). Based on the split-window algorithm proposed by Qin et al. [63], Rozenstein et al. [64] developed a two-factor split-window algorithm for calculating LST using Landsat 8 data. Jiménez-Muñoz et al. [65] also developed single-channel and split-window algorithms for Landsat 8 LST inversion. However, according to the notice issued by the United States Geological Survey (USGS) [66], the calibration of Landsat 8 TIR bands indicated that the root mean square error (RMS) of Band 11 is much greater than that of Band 10. Therefore, given the greater uncertainty associated with TIR Band 11, the USGS recommends using TIR Band 10 to retrieve LST independently. For previous Landsat TM or ETM+ data with a single TIR band, Qin et al.'s [45] mono-window algorithm is widely used for LST inversion as it only requires three additional parameters, namely the effective mean atmospheric temperature, atmospheric transmittance, and land surface emissivity. Based on this, Wang et al. [66] developed an improved mono-window (IMW) algorithm specifically for Landsat 8 data, which has been shown to be highly accurate in its application. The IMW algorithm can be expressed as follows: where T sur is the LST; a and b are the linearized regression coefficient of Planck blackbody radiation function for TIR Band 10 (for the temperature range 0-50 • C: a = −62.7182 and b = 0.4339); T 10 is the brightness temperature calculated using TIR Band 10 [66]; T air_e is the effective mean atmospheric temperature; and C and D are intermediate parameters determined using atmospheric transmittance (τ) and land surface emissivity (ε), respectively. As Xuzhou is located in the middle latitudes, T air_e can be calculated based on its linear relationship with near-surface air temperature (T air ) during the summer (Equation (46)) and winter (Equation (47)) [45]: The emissivity of a water body in the thermal band is very high (close to a blackbody) and, therefore, water emissivity can be estimated directly as ε w = 0.991 [66]. For natural (vegetation and soil) and anthropogenic (impervious surfaces) land surface types, land surface emissivity can be estimated as follows [45]: where P v is the vegetation proportion of a mixed pixel calculated from the NDVI; R v is the pixel radiation ratio of vegetation, where R i refers to R s and R m , which are the pixel radiation ratios of soil and anthropogenic land, respectively; ε v is the land surface emissivity of a completely vegetated pixel (ε v = 0.973); ε i refers to ε s and ε m , which are the emissivity of pure soil (ε s = 0.966) and anthropogenic pixel (ε m = 0.962) surfaces, respectively; and d ε is the thermal radiation interaction between vegetation and soil (or anthropogenic surfaces), which can be estimated using P v [66]. As atmospheric transmittance (τ) is mainly affected by atmospheric water content (w) [45], Wang et al. [66] developed regression equations for Landsat 8 TIR Band 10 using the Moderate Resolution Atmospheric Transmission simulation program (MODTRAN 4). For different ranges of w during the mid-latitude summer and winter, τ can thus be estimated using the equations in Table 3. As water content (w) data have not been directly observed in the study area, the empirical algorithm derived from the observation data of 28 Chinese meteorological stations was applied to estimate water content [33].
where e is the atmospheric water vapor pressure; ϕ = 34.24 • N is the central latitude of the study area; and H = 40 m is the average elevation of the study area.

Accuracy of Modeled ET
The representative spatial range of the ground flux observations was determined from the observation height, atmospheric stability, wind direction, wind speed, and surface coverage. Thus, if ground flux observations are directly used to verify the inversion values for the corresponding pixels, large errors are likely to occur [67]. Therefore, the footprint model [67,68] was used to evaluate the accuracy of the ET estimations. In the footprint model, the region that dominates the contribution rate of ground flux observations is called the flux source area (Ω ). The model describes the spatial distribution relationship between ground flux observations and the source areas at or near the Earth surface ( Figure 5), which can be expressed as follows: where is the downwind distance; is the crosswind distance; is the observation height of EC; ( , ) is the lateral concentration diffusion function; and ( , ) is the crosswind integral function [69].

Accuracy of Modeled ET
The representative spatial range of the ground flux observations was determined from the observation height, atmospheric stability, wind direction, wind speed, and surface coverage. Thus, if ground flux observations are directly used to verify the inversion values for the corresponding pixels, large errors are likely to occur [67]. Therefore, the footprint model [67,68] was used to evaluate the accuracy of the ET estimations. In the footprint model, the region that dominates the contribution rate of ground flux observations is called the flux source area (Ω p ). The model describes the spatial distribution relationship between ground flux observations and the source areas at or near the Earth surface ( Figure 5), which can be expressed as follows: where x is the downwind distance; y is the crosswind distance; z m is the observation height of EC; D y (x, y) is the lateral concentration diffusion function; and f y (x, z m ) is the crosswind integral function [69].
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 22 As shown in Equation (53) and Figure 5b, considering the area flux contribution, the weighted sum of modeled ET values ( ) of all the pixels covered by the source area were converted to source ET values ( ), which have the same spatial representativeness as the values obtained using the EC ( ). Therefore, can be directly validated using , as shown in Figure 6.   Figure 6b shows that points are evenly distributed near to the 1:1 line, and Table 4 reveals that at the confidence level of 99%, their linear regression also shows a good fit (R 2 = 0.8691) and a strong positive correlation (r = 0.9322) with P＜0.01. In addition, the mean relative error (MRE) and the root mean square error (RMSE) between and are 0.1976 and 27.15 W·m -2 , respectively. These results indicate that the accuracy of the modeled ET values is acceptable, and the urban RS-PM model is, therefore, considered reliable for simulating urban ET.  As shown in Equation (53) and Figure 5b, considering the area flux contribution, the weighted sum of modeled ET values (ET i ) of all the pixels covered by the source area were converted to source ET values (ET s ), which have the same spatial representativeness as the values obtained using the EC (ET o ). Therefore, ET s can be directly validated using ET o , as shown in Figure 6.
As shown in Equation (53) and Figure 5b, considering the area flux contribution, the weighted sum of modeled ET values ( ) of all the pixels covered by the source area were converted to source ET values ( ), which have the same spatial representativeness as the values obtained using the EC ( ). Therefore, can be directly validated using , as shown in Figure 6.   Figure 6b shows that points are evenly distributed near to the 1:1 line, and Table 4 reveals that at the confidence level of 99%, their linear regression also shows a good fit (R 2 = 0.8691) and a strong positive correlation (r = 0.9322) with P＜0.01. In addition, the mean relative error (MRE) and the root mean square error (RMSE) between and are 0.1976 and 27.15 W·m -2 , respectively. These results indicate that the accuracy of the modeled ET values is acceptable, and the urban RS-PM model is, therefore, considered reliable for simulating urban ET.    Figure 6b shows that points are evenly distributed near to the 1:1 line, and Table 4 reveals that at the confidence level of 99%, their linear regression also shows a good fit (R 2 = 0.8691) and a strong positive correlation (r = 0.9322) with P < 0.01. In addition, the mean relative error (MRE) and the root mean square error (RMSE) between ET s and ET o are 0.1976 and 27.15 W·m −2 , respectively. These results indicate that the accuracy of the modeled ET values is acceptable, and the urban RS-PM model is, therefore, considered reliable for simulating urban ET.

Relationship between Urban ET and LST
The ET and LST values of all pixels (excluding water) in the study area were extracted and plotted in the scatter plots shown in Figure 7. The relationship between these two variables was evaluated using correlation and linear regression analyses, also shown in Figure 7 and Table 5. The correlation significance test indicated that at the confidence level of 99%, the negative correlation between ET and LST is significant for each date (P < 0.01). Moreover, urban ET and LST were more strongly correlated between May and September (|r| > 0.5 and R 2 > 0.3) but were more weakly correlated between November and March (|r| < 0.5 and R 2 < 0.3). It should be noted that the significance of the correlation between ET and LST in October depends on the weather conditions in any single year. This indicates that urban ET had a stronger impact on LST during the warm months and a weaker influence during the cold months.

Relationship between Urban ET and LST
The ET and LST values of all pixels (excluding water) in the study area were extracted and plotted in the scatter plots shown in Figure 7. The relationship between these two variables was evaluated using correlation and linear regression analyses, also shown in Figure 7 and Table 5. The correlation significance test indicated that at the confidence level of 99%, the negative correlation between ET and LST is significant for each date (P＜0.01). Moreover, urban ET and LST were more strongly correlated between May and September (|r| > 0.5 and R 2 > 0.3) but were more weakly correlated between November and March (|r| < 0.5 and R 2 < 0.3). It should be noted that the significance of the correlation between ET and LST in October depends on the weather conditions in any single year. This indicates that urban ET had a stronger impact on LST during the warm months and a weaker influence during the cold months.

Relationship between the Intensity of Urban ET and the UHI Effect
To further explore the impact of urban ET on the UHI, the ET and LST values for each period were normalized (Equations (54) and (55)) using ET intensity (ETi) and UHI intensity (UHIi). ETi and UHIi were classified into five levels based on the following conditions: (1) All pixels were sorted in ascending order according to their intensity value (ETi or UHIi, respectively), and (2) the first 20% of pixels were categorized as level 1 (LV 1), and then each subsequent 20% of pixels were categorized as level 2 to level 5 (LV 2 to LV 5), respectively. Finally, the different levels of ETi and UHIi were visualized and, by overlaying them, the proportion of each UHIi level at each ETi level was assessed (Figure 8).

Relationship between the Intensity of Urban ET and the UHI Effect
To further explore the impact of urban ET on the UHI, the ET and LST values for each period were normalized (Equations (54) and (55)) using ET intensity (ET i ) and UHI intensity (UHI i ). ET i and UHI i were classified into five levels based on the following conditions: (1) All pixels were sorted in ascending order according to their intensity value (ET i or UHI i , respectively), and (2) the first 20% of pixels were categorized as level 1 (LV 1), and then each subsequent 20% of pixels were categorized as level 2 to level 5 (LV 2 to LV 5), respectively. Finally, the different levels of ET i and UHI i were visualized and, by overlaying them, the proportion of each UHI i level at each ET i level was assessed (Figure 8).
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing May and October was significantly higher than between November and March. Finally, the sum of the areas classified as UHIi = LV 1 and LV 2 are highest in those regions where ETi = LV 5 (61%-85%, except for 2015-12-21 and 2016-11-05). This indicates that the highest ET intensities had the most significant impact on the UHI. In other words, ET played an important role in regulating the UHI in those regions where the intensity of ET was highest during the summer (or near summer) period.

The Effect of High ET Intensity on UHI
Given our finding that higher ET intensities had the strongest influence on the UHI effect, especially during the warmer months of the year, it is necessary to clarify the mechanism underpinning this relationship. To do this, we used buffer analysis. First, data from six periods between May and October were selected, and those regions where ETi = LV 5 were extracted. Second, using the LV 5 areas as the central areas, five buffer zones were established at 30 m intervals (the distance of one pixel), as shown in Figure 9. Finally, the vector layers of the core areas (ETi = LV 5) and the five buffer zones were overlaid onto the corresponding ET and LST images, and the average ET and LST values for each layer were calculated ( Figure 10). It should be noted that to avoid the interference of waterbodies, their ET and LST values were assigned as null values so that they were excluded from the analysis. Three important features can be noted in Figure 8. First, with an increase in the ET intensity level, the UHI intensity level generally decreased, which indicates that the UHI intensity is weaker in regions with high ET intensity levels. Second, the impact of ET intensity on the UHI intensity between May and October was significantly higher than between November and March. Finally, the sum of the areas classified as UHI i = LV 1 and LV 2 are highest in those regions where ET i = LV 5 (61%-85%, except for 2015-12-21 and 2016-11-05). This indicates that the highest ET intensities had the most significant impact on the UHI. In other words, ET played an important role in regulating the UHI in those regions where the intensity of ET was highest during the summer (or near summer) period.
where T sur,max and ET max are the maximum values of LST and ET, respectively; T sur,min and ET min are the minimum values of LST and ET, respectively.

The Effect of High ET Intensity on UHI
Given our finding that higher ET intensities had the strongest influence on the UHI effect, especially during the warmer months of the year, it is necessary to clarify the mechanism underpinning this relationship. To do this, we used buffer analysis. First, data from six periods between May and October were selected, and those regions where ET i = LV 5 were extracted. Second, using the LV 5 areas as the central areas, five buffer zones were established at 30 m intervals (the distance of one pixel), as shown in Figure 9. Finally, the vector layers of the core areas (ET i = LV 5) and the five buffer zones were overlaid onto the corresponding ET and LST images, and the average ET and LST values for each layer were calculated ( Figure 10). It should be noted that to avoid the interference of waterbodies, their ET and LST values were assigned as null values so that they were excluded from the analysis.    Figure 10 shows that the average ET values in the core zones are significantly higher than in each successive (outward) buffer zone, while the average LST values show the opposite trend. Indeed, the average surface ET values decrease outwards from the core zones while the average LST values increase. Furthermore, areas further away from the core zones show lower rates of change in ET (decreases) and LST (increases). These characteristics indicate that regions with high ET intensities exert a regulatory effect on the thermal environment of their surrounding area and that this relationship weakens with distance. To further quantify the regulating effect of areas with high ET intensities on their surrounding areas, the variation of ET (∆ ) and LST (∆ ) between the adjacent buffer layers was evaluated, as shown in Table 6.      Figure 10 shows that the average ET values in the core zones are significantly higher than in each successive (outward) buffer zone, while the average LST values show the opposite trend. Indeed, the average surface ET values decrease outwards from the core zones while the average LST values increase. Furthermore, areas further away from the core zones show lower rates of change in ET (decreases) and LST (increases). These characteristics indicate that regions with high ET intensities exert a regulatory effect on the thermal environment of their surrounding area and that this relationship weakens with distance. To further quantify the regulating effect of areas with high ET intensities on their surrounding areas, the variation of ET (∆ ) and LST (∆ ) between the adjacent buffer layers was evaluated, as shown in Table 6.   Figure 10 shows that the average ET values in the core zones are significantly higher than in each successive (outward) buffer zone, while the average LST values show the opposite trend. Indeed, the average surface ET values decrease outwards from the core zones while the average LST values increase. Furthermore, areas further away from the core zones show lower rates of change in ET (decreases) and LST (increases). These characteristics indicate that regions with high ET intensities exert a regulatory effect on the thermal environment of their surrounding area and that this relationship weakens with distance. To further quantify the regulating effect of areas with high ET intensities on their surrounding areas, the variation of ET (∆ET Abu f ) and LST (∆T Abu f ) between the adjacent buffer layers was evaluated, as shown in Table 6.
From the outermost buffer zone to the innermost buffer zone, higher increases in ∆ET Abu f corresponded to greater decreases in ∆T Abu f . Normalized ∆ET Abu f and ∆T Abu f are also shown in Figure 11a, which indicates that variations in these two parameters are highly coincident. At the confidence level of 99%, linear regression analysis between ∆ET Abu f and ∆T Abu f (Figure 11b and Table 7) gave a Pearson correlation coefficient (r) of −0.9628 and an R 2 of 0.9270 with P < 0.01, indicating a strong negative linear correlation. Our analyses, therefore, indicate that those urban regions with high ET intensities play an important role in alleviating temperatures in the surrounding environment. According to the linear fitting equation, in high ET intensity areas, for every 10 W·m −2 increase in ET intensity, LST decreases by approximately 0.56 K within a 150 m radius (excluding other cooling factors). From the outermost buffer zone to the innermost buffer zone, higher increases in ∆ corresponded to greater decreases in ∆ . Normalized ∆ and ∆ are also shown in Figure 11a, which indicates that variations in these two parameters are highly coincident. At the confidence level of 99%, linear regression analysis between ∆ and ∆ (Figure 11b and Table 7) gave a Pearson correlation coefficient (r) of -0.9628 and an R 2 of 0.9270 with P＜0.01, indicating a strong negative linear correlation. Our analyses, therefore, indicate that those urban regions with high ET intensities play an important role in alleviating temperatures in the surrounding environment. According to the linear fitting equation, in high ET intensity areas, for every 10 W·m −2 increase in ET intensity, LST decreases by approximately 0.56 K within a 150 m radius (excluding other cooling factors).

Discussion
Our verification results for the urban RS-PM model [35] show a good inversion accuracy was achieved and that the parameters of the model are easy to obtain. Compared with the modified multisource parallel model [33], which can only be applied to remote sensing data with multiple TIR bands, the urban RS-PM model has better applicability and generalizability for urban ET estimation.
Most previous studies of mitigating the UHI effect have focused on factors including vegetation indices [19,71,72], the "cool island effect" of urban parks (which can be affected by urban parks' size, land surface coverage types and shape, etc.) [73][74][75], and land surface moisture [76]. By approaching the challenge from a new perspective, we explored the magnitude and spatial characteristics of the

Discussion
Our verification results for the urban RS-PM model [35] show a good inversion accuracy was achieved and that the parameters of the model are easy to obtain. Compared with the modified multi-source parallel model [33], which can only be applied to remote sensing data with multiple TIR bands, the urban RS-PM model has better applicability and generalizability for urban ET estimation.
Most previous studies of mitigating the UHI effect have focused on factors including vegetation indices [19,71,72], the "cool island effect" of urban parks (which can be affected by urban parks' size, land surface coverage types and shape, etc.) [73][74][75], and land surface moisture [76]. By approaching the challenge from a new perspective, we explored the magnitude and spatial characteristics of the urban ET alleviation effect on the UHI.
Our results indicate that ET from urban vegetation and bare soil has a mitigating effect on the UHI; however, this effect is not always significant. During the colder months of the year (November-March) no significant correlation was found between ET and LST (|r| = 0.1058-0.4863 and R 2 = 0.0111-0.2334) (Figure 7), Furthermore, during these months, increases in the intensity of ET were not associated with significant changes in the intensity of the UHI (Figure 8). For example, for 5 November 2016 and 21 December 2015, the proportions of areas with a low UHI intensity (LV 1 and LV 2) occurring in areas with the highest ET intensities (LV 5) were approximately 47.4% and 47.3%, respectively, which is notably lower than the warmer months of the year. These differences are explained by the fact that ET rates are much lower in winter than in summer (Figure 4), meaning that urban ET is less important in regulating the UHI during the cooler autumn/winter period.
Compared with the colder months of the year, the alleviating effect of urban ET on the heat island effect was found to be highly significant during the warmest months; however, this effect was strongly affected by the intensity of ET. As shown in Figure 8, for each category of decrease in ET intensity between May and October, the area corresponding to a weak UHI effect significantly decreased. Our findings suggest that when the intensity of ET decreased from LV 5 to LV 4, the average rate of decrease in the area experiencing a LV 1 UHI intensity (the lowest intensity level) was 59.9%. Therefore, only high ET intensities were effective at alleviating the UHI effect.
We have also shown that the regulating influence of high-intensity ET areas on the UHI decreases rapidly with distance ( Figure 10). Our correlation and fitting results ( Figure 11) indicate that the degree of variation in ET is a dominant factor leading to this phenomenon. As the major driver of urban ET is vegetation, factors such as vegetation shading [77,78], size [13,79], and type [75] are known to be important controls on the cooling effect. Studies on the cooling effects of urban parks [10,80] have also found that this effect is limited to distances of 35 to 840 m from parks [81][82][83]. Therefore, the magnitude of the thermal regulating effect of high-intensity ET areas is dependent on the ET intensity alongside the characteristics of the vegetation.

Conclusions
Based on the modified RS-PM model, ET in urban areas of Xuzhou was successfully inversed for the period 2014-2018 using Landsat 8 data. The relationship between urban ET and LST was significantly negative during the warmer months of the year (May-October), but during the colder months (November-March), there was no significant correlation. This was further verified through the spatial superposition analysis of the intensity of ET and UHI, demonstrating that high-intensity regions had the most significant alleviating effect on the UHI. According to our buffer analysis, the variation in ET is a key factor regulating the UHI during the warmer summer months. We found that for every 10 W·m −2 increase in ET, LST decreased by 0.56 K. Our findings provide a new perspective for the improvement of urban thermal comfort, which can be applied to urban management, planning, and natural design.
To explore these relationships between ET and UHI in more detail, further research is required on the cooling effects of ET with different vegetation characteristics (such as shading capacity, density of planting, and size), and the application of high-resolution remote sensing data could be especially advantageous.
Author Contributions: Conceptualization, Y.W. and Y.Z.; methodology, K.Q. and Y.Z.; software, Y.Z.; validation, N.D. and K.Q.; formal analysis, N.D. and X.Y.; writing-original draft preparation, Y.W.; writing-review and editing, X.Y.; funding acquisition, N.D. and X.Y. All authors have read and agreed to the published version of the manuscript.