Combined Effects of Impervious Surface Change and Large-Scale Afforestation on the Surface Urban Heat Island Intensity of Beijing, China Based on Remote Sensing Analysis

: Urban heat island (UHI) attenuation is an essential aspect for maintaining environmental sustainability at a local, regional, and global scale. Although impervious surfaces (IS) and green spaces have been confirmed to have a dominant effect on the spatial differentiation of the urban land surface temperature (LST), comprehensive temporal and quantitative analysis of their combined effects on LST and surface urban heat island intensity (SUHII) changes is still partly lacking. This study took the plain area of Beijing, China as an example. Here, rapid urbanization and a large-scale afforestation project have caused distinct IS and vegetation cover changes within a small range of years. Based on 8 scenes of Landsat 5 TM/7ETM/8OLI images (30 m × 30 m spatial resolution), 920 scenes of EOS-Aqua-MODIS LST images (1 km × 1 km spatial resolution), and other data/information collected by different approaches, this study characterized the interrelationship of the impervious surface area (ISA) dynamic, forest cover increase, and LST and SUHII changes in Beijing’s plain area during 2009–2018. An innovative controlled regression analysis and scenario prediction method was used to identify the contribution of ISA change and afforestation to SUHII changes. The results showed that percent ISA and forest cover increased by 6.6 and 10.0, respectively, during 2009–2018. SUHIIs had significant rising tendencies during the decade, according to the time division of warm season days (summer days included) and cold season nights (winter nights included). LST changes during warm season days responded positively to a regionalized ISA increase and negatively to a regionalized forest cover increase. However, during cold season nights, LST changes responded negatively to a slight regionalized ISA increase, but positively to an extensive regionalized ISA increase, and LST variations responded negatively to a regionalized forest cover increase. The effect of vegetation cooling was weaker than ISA warming on warm season days, but the effect of vegetation cooling was similar to that of ISA during cold season nights. When it was assumed that LST variations were only caused by the combined effects of ISA changes and the planting project, it was found that 82.9% of the SUHII rise on warm season days (and 73.6% on summer days) was induced by the planting project, while 80.6% of the SUHII increase during cold season nights (and 78.9% during winter nights) was caused by ISA change. The study presents novel insights on UHI alleviation concerning IS and green space planning, e.g., the importance of the joint planning of IS and green spaces, season-oriented UHI mitigation, and considering the thresholds of regional IS expansion in relation to LST changes.


Introduction
The urban heat island (UHI) effect occurs when urban areas have higher temperatures than surrounding rural areas [1], which can contribute to a series of environmental and socio-economic problems, e.g., increasing the frequency and persistence of heat waves [2,3], raising energy consumption [4], accelerating air pollution [5], and threatening public health [6]. In the current context of global climate change and rapid urbanization, UHI attenuation based on a better understanding of its influential factors and mechanisms is a pressing issue for sustainable development [7].
UHI can be measured via both the atmospheric and land surface temperature. The latter, in terms of surface UHI (SUHI), has been studied more frequently due to the advantage of data acquisition through remote sensing imageries [8]. Urban land surface characteristics (e.g., albedo, emissivity, and thermal capacity) that are mostly determined by land-use/land-cover (LULC) varieties have been found to have a major impact on the SUHI [9]. Accordingly, a number of studies have focused on discovering the contributions of both the LULC composition and configuration to the land surface temperature (LST). Commonly, on hot summer days, impervious surfaces (IS) increase LST by raising sensible heat flux, whereas green spaces (especially forests) lower the ambient air temperature through reducing the absorption of solar radiance and increasing latent heat flux [10,11]. Water bodies also have negative effects on LST due to their high heat capacity and evaporative cooling [12]. Additionally, the effects of the LULC configuration on LST variations have been confirmed, e.g., clustered impervious surfaces/vegetation cover act more strongly in elevating/lowering LST [13], and the rational use of connectivity from a landscape source-sink point of view can bring additional cooling effects [14]. On the whole, the composition of LULC has been proven to be more important in influencing LST than the configuration based on macro-scale remote sensing observations, and the proportions of impervious surface (IS) and green spaces have been shown to have a dominant effect on LST spatial differentiation in different seasons [15][16][17][18][19][20][21].
Compared with the pattern analysis of LST in relation to the LULC composition and configuration, the temporal dynamic of LST in response to LULC conversions seems to be more complicated and has been given less attention. Yu et al. [22] explored summertime LST changes in an urban agglomeration of southern China. Their study concluded that the transformation of built-up areas from all other land cover types can increase LST, except for bare land, and the transformation of forestland from all other land cover types decreases LST. Akinvemi et al. [23] and Muro et al. [24] studied land cover variations of dryland and wetland regions in relation to LST changes, respectively, also concluding that positive LST trends corresponded to deforestation and farmland expansion, while negative LST trends corresponded to forestation processes. In addition to studies following the land use transfer matrix, research has looked into the relationships between individual LULC component (especially IS and vegetation cover) dynamics and LST changes. For example, time-series models were used to measure and predict LST changes in relation to urban expansion [25][26][27][28]. It is anticipated that urban expansion-induced global warming between 2015 and 2050 will be about half to two times as strong as that caused by greenhouse gas emissions [29]. Vegetation cooling has been valued as a nature-based solution (NBS) for counteracting the negative effects of urbanization [30][31][32]. However, the effects of IS and vegetation cover changes of the same area in influencing LST are inequivalent, as IS change generally has a stronger relationship with LST [33][34][35]. Moreover, the cooling/warming effect and magnitude of vegetation gain/loss are still debated, as few studies have found that they highly depend on the growing stages, previous land cover types, and biophysical characteristics [22,[36][37][38].
The SUHI intensity (SUHII) reflects the UHI effect directly and quantitatively. So far, three types of methods have been applied to estimate SUHII: LST as a proxy of SUHI; the LST difference between the urban and reference areas; and non-parametric models [39,40]. While the second method is the most commonly used, the third of non-parametric models can avoid the problem of defining urbanrural boundaries that cause additional SUHII deviations [41][42][43]. Research on the temporal dynamic analysis of SUHII in response to LULC changes has been very limited. Li et al. [44] studied approximately 5000 urban areas of the conterminous United States and showed that the urban area size alone can explain as much as 87% of the variance of SUHII during urban expansion. Yao et al. [45] found that rural greening contributed to 22.5% of the increased daytime SUHII in the growing season at a global scale during 2001-2017, indicating that afforestation or reforestation could also exacerbate the UHI effect without rational spatial planning. Overall, despite the high-ranking importance of IS and vegetation cover changes in influencing urban LST having been demonstrated, their effects on SUHII changes remain largely unknown.
In general, past research has provided plentiful references for urban landscape planning in terms of coping with the UHI effect, especially IS and green spaces at a local scale. As shown by Zhou et al. [39,40], SUHII publications have grown exponentially since 2005 due to the growing interest associated with rapid urbanization and global warming, in addition to the advancement of new remote sensing techniques. However, with respect to these pieces of literature, debate exists, for example, on aspects of the cooling effects of afforestation or reforestation, as well as their role in SUHII attenuation. Moreover, although the changing trend of LST in response to LULC conversions has been frequently discussed, quantitative relationships are still unclear due to inadequate timeseries analysis. Furthermore, many relevant studies have focused on a specific time (generally summer), whereas inter-annual variations of LST changes in response to LULC conversions are yet to be explored in seasonal and diurnal divisions. Additionally, focus on practically implemented urban strategies could refresh the conclusions of this field of study, especially if advanced remote sensing techniques were applied [46]. It is therefore important to follow opportunities that could motivate new insights on this "old" topic.
Due to rapid development, some Chinese megacities (e.g., Beijing) have seen extensive IS expansion and greening efforts during the past two or three decades [47]. From 2012 to 2015, the Beijing municipal government accomplished a large-scale planting project to increase forest cover by over 10% in its plain area-the One-Million-Mu (666 km 2 ) Plain Afforestation Project [48]. This unprecedented greening project motivated the conception of this study, leading to three specific questions: (a) How have IS, forest cover, and LST changed over the range of years covering the largescale afforestation period in the plain area of Beijing?; (b) what is the interrelationship among IS, forest cover, and LST changes with seasonal and diurnal variations?; and (c) how have IS change and afforestation influenced the SUHIIs of this area? Based on remote sensing data of Landsat 5 TM/7ETM/8OLI images (30 m × 30 m spatial resolution), EOS-Aqua-MODIS LST images (1 km × 1 km spatial resolution), and other data/information collected from different approaches, a series of methods, including temporal trend analysis, regression analysis, and a new technique of scenario prediction, were used to find answers to the three questions. The study aimed to provide implications for sustainable urban landscape planning in coping with UHI, based on a novel understanding of the combined effects of IS and forest cover changes on LST and SUHII variations.

Study Area
Beijing covers an area of 16,410 km 2 . It has a warm temperate and monsoon-influenced continental climate, with hot humid summers and cold, windy, and dry winters. The Taihang and Yanshan Mountains shield the city to the north, northwest, and west, while the North China Plain opens to the south and east of the city (Figure 1b). The plain area (6338 km 2 , altitude generally lower than 100 m) hosts the majority of the population of Beijing and has experienced drastic urban expansion driven by concentric development, with the appearance of five city ring roads ( Figure 1c) [49,50]. A negative influence of urban expansion on the thermal environment of Beijing has been found [51,52]. From 1960 to 2008, the annual mean temperature of Beijing increased 0.39 °C per year on average, which is much higher than the global warming rate of 0.13 °C [53]. A previous study also revealed that latent heat flux (LHF) in the plain areas changed more drastically from 1997 to 2017 compared with the mountain areas [51].  Image preprocessing was conducted with the ENVI software (version 5.2; Exelis VIS, Inc., Boulder, CO, USA). Procedures included radiometric correction, FLASH atmospheric correction, image mosaicing, and image cutting. The Modified Normalized Difference Water Index (MNDWI) was calculated following Equation (1) [54], in order to extract water cover by exploring and setting threshold values.

Data Collection and Preprocessing
where GREEN and SWIR represent the reflectance of green (band 2 of Landsat 5 TM and 7 ETM+, and band 3 of Landsat 8 OLI/TIRS) and shortwave infrared (SWIR 1, band 5 of Landsat 5 TM and 7 ETM+, and band 6 of Landsat 8 OLI/TIRS) spectral bands, respectively.
Furthermore, the Biophysical Composition Index (BCI) was calculated to extract ISA by exploring and setting threshold values. BCI was proposed by Deng and Wu [55] following Ridd's [56] conceptual vegetation-impervious surface-soil triangle model, and was derived by reexamining the components of Tasseled Cap (TC) transformation using Equation (2): where H is "high albedo", the normalized TC1; L is "low albedo", the normalized TC3; and V is "vegetation", the normalized TC2. These three factors can be calculated as follows: where TCi (i = 1, 2, and 3) represents the first three TC components, and TCimin and TCimax are the minimum and maximum values of the ith TC component, respectively. TC transformation of Landsat 5 TM and Landsat 7 ETM+ images was realized using the default parameters in ENVI 5.2, while TC transformation of Landsat 8 OLI/TIRS images was conducted following the study of Baig et al. [57]. Finally, 300 checking points were randomly selected to testify the classification results in reference to high-resolution Google Earth images.

Calculation of the Diurnal and Seasonal Average LST and SUHII
MODIS LST data have been frequently used in UHI studies due to their high temporal resolution [58]. In this study, time series LST data were obtained from version 6 EOS-Aqua-MODIS LST data (MOD11A2, 8-day composite, 1 km × 1 km spatial resolution) through the Google Earth Engine [59]. We chose the Aqua satellite as it passes over Beijing at local time of 01:30 (nighttime) and 13:30 (daytime), which is closer to the lowest and highest temperature during the day compared with the Terra satellite at the local time A non-parametric model of the linear regression function between LST and regionalized ISA to calculate SUHII proposed by Li et al. [60] was used in this study ( Figure 2). The steps applied were as follows. Firstly, a 1 km × 1 km fishnet was created within the study area considering the resolution of MODIS LST data. The proportion of water cover and ISA in each grid cell of the fishnet in the years of 2009, 2011, 2014, and 2017 was calculated. Secondly, grid cells with more than half of their total area outside the boundary and those with a water coverage higher than 0.25 were excluded. This meant that 6325 grid cells were left. Thirdly, grid cells were converted into points using the feature to point module in the Arcgis software (version 10.5; Esri, Inc., RedLands, CA, USA), and the ISA proportion was regionalized using the Kernel Density Estimation (KDE) method. The kernel radius was 3 km and the resolution of the output raster cell was 1 km. Fourth, all the output pixels within the study area were clustered based on 1% intervals of ISAKDE. The values of ISAKDE and LST of the same cluster were averaged to counteract the sub-interval variation of LST caused by the heterogeneity of the urban surface. Finally, the linear regression slopes of percentile-averaged ISAKDE and LST in the years from 2009 to 2018 were calculated and defined as SUHIIs if the coefficient of determination (R 2 ) employed to evaluate the accuracy of the function was higher than 0.   [60]. LSTr means land surface temperature (LST) in the rural areas, and ISAKDE means the regionalized impervious surface area (ISA) obtained using the Kernel Density Estimation (KDE) method.

Temporal Trend Analysis of SUHII and Pixel-Wise LST
The yearly value of SUHIIs and pixel-wise LST was averaged in three periods, i.e., Period I (2009-2011, before the planting), Period II (2012-2015, during the planting), and Period III (2016-2018, after the planting), in order to find potential rules for the changes related to the planting. Simultaneously, the temporal trends of SUHII and pixel-wise LST in time series of 2009-2018 were detected by using the Mann-Kendall (MK) [61,62] and Sen's slope estimator [63] non-parametric tests to further detect the rules. The MK test was used to disclose the statistical significance of the trend in time series, and the Sen's slope estimation was employed to discover the trend magnitude. These two methods have been widely used in climate and environmental research because no distributional hypotheses for the variables are required [28,64]. The calculation of the MK and Sen's slope test was conducted using the mkttest function of the modifiedmk package, and the sqmk function of trendchange package with R (version 3.3.1; R Development Core Team, Vienna, Austria). The visualization of pixel-wise variables was realized in Arcgis 10.5.

Regression Analysis of LST Changes in Response to the ISA Dynamic and Afforestation
Basic data on the variables, including the LST, ISAKDE, and new plantation areas in the 6325 grid cells (1 km*1 km) covering the study area described above in Section 2.2.3 were collected. The average LST of Period I and III in each grid and its increase from Period I to III were calculated, and recorded as ΔLST. Similarly, the increase of ISAKDE from Period I (ISAKDE2011) to Period III (ISAKDE2017) in each grid was computed, and recorded as ΔISAKDE. The area of new plantations in each grid was also regionalized using the Kernel Density Estimation (KDE) method following the same steps in Section 2.2.3 and was named FORKDE.
The ordinary least squares (OLS) model was used as a conventional regression model to estimate geographical issues [65]. It was adopted in this study to figure out the general effects of ISA change and large-scale afforestation on LST changes from a global view of the study area. The process of OLS estimation was realized using the Ordinary Least Squares module in Arcgis 10.5. Furthermore, in order to unveil the individual contribution of ISA change and afforestation to SUHII change, the data were divided into two sub-datasets of "ISA unchanged" and "ISA changed", as well as "without new plantations" and "with new plantations", following the steps displayed in the supplementary file ( Figure S1). To counteract the sub-interval variations caused by spatial heterogeneity, processes to obtain the pixel group and data average were conducted following the 1% intervals clustering method, as described in Section 2.2.3. A linear regression simulation was chosen as the preferential model and was realized using Microsoft Excel 2016. For scatter diagrams that showed evident trend transition and inflection points, a piecewise linear regression simulation was conducted to improve the model reliability using the segmented function of the segmented package with R software.

Scenario Prediction
As previous studies have proved that artificial surfaces and green space determine LST variations in all seasons [17,20], it was hypothesized that LST change was only caused by the combined effects of ISA change and the planting project. To unveil the share of contributions of ISA change and afforestation to SUHII variation in different time divisions, the controlled regression models mentioned in Section 2.3.2 using sub-datasets of "ISA unchanged" and "without new plantations" were employed. Regression models of higher correlation coefficients were selected to predict ΔLST from Period I to Period III under two scenarios: (A) Without the afforestation project, and (B) without ISA change during the transition from Period I to Period III. The actual SUHII of different time divisions in Period I and Period III, as well as the predicted SUHII in Period III under the two scenarios, was calculated. Finally, the actual and predicted SUHII increases from Period I to Period III were compared, and the difference could be attributed to the assumption of "without the afforestation project" or "without ISA change". Therefore, the contribution shares were uncovered.

ISA Dynamics, Location of the Planting Sites, and SUHII Changes
The accuracy of land cover classification of 2009, 2011, 2014, and 2017 was 87.9%, 85.3%, 88.2%, and 82.3%, respectively. The ISA percentage in the plain area of Beijing consistently increased from 38.1% in 2009 to 44.7% in 2017 ( Figure 3). The ISA increment seems to be more intense near and outside the Sixth Ring Road. When urban core areas and suburb/rural areas were zoned by the contour line of ISAKDE of 50% [28], the ratio of urban core areas increased from 1752 km 2 in 2009 to 2075 km 2 in 2017. The expansion of urban core areas seems to be more obvious in northeastward and southward directions.   (Tables S1 and S2). The regression functions seem well-fitted in all of the time divisions, except for spring days, winter days, and cold season days, which have R 2 values generally lower than 0.8. The highest SUHII across a year occurred on summer days, ranging from 6.17 to 9.26, followed by winter nights, ranging from 5.48 to 8.00, while the lowest SUHII occurred on autumn days, varying from 2.19 to 3.96. Table 1 reveals the SUHII with a yearly rising tendency in all diurnal and seasonal divisions during 2009-2018, except for autumn days, because all of the values of ZMK are positive, in addition to autumn days. SUHII for summer days, winter nights, warm season days, and cold season nights increased significantly. Sen's slope values of SUHII for these four time divisions were higher than others, indicating higher growth rates of SUHII. Average SUHIIs in the three periods divided by the planting project showed a consistent rise in all of the time divisions (values available), except for autumn days, summer nights, and warm season nights ( Figure 5). Considering the results above, the focus in the subsequent study was on summer days, winter nights, warm season days, and cold season nights, as SUHII and its increase in these four time divisions were more evident.  (Figure 6c,e). However, during winter nights and cold season nights, the distribution was different, i.e., high-value areas are mainly inside or around the Sixth Ring Road, especially in the south, where a few clustered pixels display a significant LST increase; lowvalue areas are generally located in the northeast and southeast, and very few pixels show a significant LST decrease (Figure 6d,f). This disparity suggests different driving mechanisms of LST changes during warm season days and cold season nights.

LST Changes in Response to ISA Dynamics and Afforestation
Distributions of the average LST in the three periods are displayed in the Supplementary Materials, Figure S2. Table 2 reveals that FORKDE was significantly negatively correlated with ΔLST during summer days, winter nights, and cold season nights, but insignificantly positively correlated with ΔLST on warm season days. ΔISAKDE was significantly positively correlated with ΔLST on summer days and warm season days, and significantly negatively correlated with ΔLST during winter nights and cold season nights. On summer days, the regression coefficient between ΔISAKDE and ΔLST was nearly 20 times the coefficient between FORKDE and ΔLST, indicating a substantial difference of contribution magnitudes to ΔLST. However, during winter nights and cold season nights, the regression coefficient between ΔISAKDE and ΔLST was slightly larger than that between FORKDE and ΔLST, suggesting a similar magnitude of contribution to LST. In total, 499 pixels with a zero value of FORKDE and 468 pixels with a zero value of ΔISAKDE were found among all of the 6325 pixels (Figure 7a or Figure 8a). In the groups with zero-valued ΔISAKDE, the average ΔLST showed a downward trend as FORKDE increased in all four time divisions, which indicates that afforestation can restrict the LST increase on both warm season days and during cold season nights (Figure 7c-f). In groups with nonzero-valued ΔISAKDE, the average ΔLST also showed a downward trend as FORKDE increased on summer days and warm season days; however, during winter nights and cold season nights, there was an inflection point, after which the trend reversed. This may be caused by the difference in the positive and negative ΔISAKDE value before and after the inflection point, reflected in Figure 7b, as the LST change in the group of nonzero-valued ΔISAKDE was related to the combined effects of afforestation and ISA change. It was also noticed that on summer days and warm season days, the group with nonzero-valued ΔISAKDE generally had a higher ΔLST compared with zero-valued ΔISAKDE, while during winter nights and cold season nights, the group with nonzero-valued ΔISAKDE had lower ΔLST. This phenomenon manifests the heating effect of impervious surfaces on warm season days and the cooling effect during cold season nights, which conforms to the above OLS estimation results. Likewise, in groups with zero-valued FORKDE, the average ΔLST showed an upward trend as ΔISAKDE increased on summer days and warm season days (Figure 8c,e), and there was a turning point (ΔISAKDE equals 0.217), after which the slope increased. During winter nights and cold season nights, the average ΔLST firstly showed a downward trend along the gradient of ΔISAKDE, and then exhibited a slightly upward trend when ΔISAKDE was higher than 0.096 during winter nights and 0.034 during cold season nights (Figure 8d,f). This indicates that ISA expansion can accelerate the LST increase on warm season days, especially when ΔISAKDE exceeds 0.217, and can restrict the LST increase during cold season nights within certain thresholds. ΔLST in the groups with nonzerovalued FORKDE displayed a similar trend to the groups with zero-valued FORKDE on summer days and warm season days. During winter nights and cold season nights, however, the trends were different. The average ΔLST firstly showed an upward trend within a ΔISAKDE value of −9.5 and −9.6 during winter nights and cold season nights, respectively. Then, the trend reversed beyond the threshold value. It could be noticed that this trend transition was closely related to the average FORKDE dynamic along the ΔISAKDE gradient revealed in Figure 8b. The values of ΔLST in the group with nonzerovalued FORKDE were generally lower than in the group with zero-valued FORKDE, from which the cooling effect of afforestation during cold season nights can be derived.

Contributions of ISA Change and Afforestation to SUHII Increases
SHUIIs of the four time divisions in Period I and III, as well as predicted scenarios of Period III, were calculated using the linear regression model (Figure 9). Models selected from Section 3.3 to predict ΔLST in the two scenarios include the piecewise regression models of ΔLST and ΔISAKDE on summer days and warm season days in the groups with zero-valued FORKDE, as well as linear regression models of ΔLST and FORKDE during winter nights and cold season nights in the groups with zero-valued ΔISAKDE.
As shown in Figure 9, SUHII on summer days increased by 1.40 from Period I (7.47) to III (8.87), whereas, without the planting project, the increase should have been 0.37. This indicates that the ISA change contributed 26.4% of the total SUHII increase on summer days, and the remaining 73.6% increase was caused by afforestation. Similarly, on warm season days, the ISA change contributed 17.1% of the SUHII increase, and the remaining 82.9% increase was caused by afforestation (Table 3). SUHII during winter nights increased by 0.98 from Period I (6.40) to III (7.38), whereas, without ISA change, the increase should have been 0. 19. This suggests that afforestation led to 19.4% of the SUHII increase during winter nights, and the remaining 80.6% increase was induced by ISA change. Similarly, during cold season nights, 19.4% of the SUHII increase was caused by afforestation and the remaining 78.9% was attributed to ISA change. Altogether, afforestation had the dominant effect on increasing SUHII on warm season days, while ISA change had the major effect on increasing SUHII during cold season nights. Figure 9. Calculation of SUHIIs in Period I, Period III, and scenarios predicted in Period III: (a) SUHII calculation during summer days; (b) SUHII calculation during winter nights; (c) SUHII calculation during warm season days; and (d) SUHII calculation during cold season nights. Note: Scenario A represents "without the afforestation project" and is applied during summer days and warm season days, and Scenario B represents "without ISA change" and is applied during winter nights and cold season nights. Table 3. Contributions of the ISA change and afforestation to SUHII increases derived from the difference of the actual and predicted SUHII in the two scenarios from Period I to Period III.

Scenarios
Time Divisions

Spatial-Temporal Changes of ISA, Urban Forests, and SUHII in Beijing's Plain Area
The study revealed that the percent of ISA in the plain area of Beijing increased by 6.6 from 2009 to 2017. Previous studies depicting the urbanization dynamics of Beijing showed divergent results because of disparities of data sources, geographic extents, urban area definitions, and time ranges, etc. [66][67][68]. A shared phenomenon was that ISA and urban areas in Beijing had a consistent increasing tendency during the past decades, and the urban-rural fringe generally extended to the northeast and south, especially near and outside the Sixth Ring Road.
The ratio of forest cover in Beijing's plain area increased by 10.0 due to the large-scale afforestation project from 2012 to 2015. However, according to the study of Huang et al. [69], only 52.34% of the total afforested/reforested area could be recognized by Landsat remote sensors in 2015. One of the reasons for this may be that most of the new forestland was converted from cropland. These two land cover types are easily confused because of their similar reflectance. Another important reason may be that there is a growing stage in which the trees are "green" enough to be captured through remote sensing. Similarly, ecological services of the new plantations, such as environment cooling, could also have a development period, and that is why, in the present study, the ten year range was divided into three periods according to the planting stage.
The regression analysis method applied in this study to calculate SUHII avoided the problem of defining urban-rural boundaries that can cause variations in the results [41,42,70]. It is based on a strong linear relationship between LST and the percent of ISA in all the year-round seasons [60,71]. However, in this study, the correlation coefficient between LST and ISAKDE for cold season days was relatively low, which requires further discussion. These time divisions were omitted in the next steps. Seasonal and diurnal SUHIIs were characterized by higher values on summer days and during winter nights based on this study, which is in line with previous studies [41,72]. Moreover, it was discovered that during warm season days and cold season nights, SHUII in Beijing's plain area showed a significant rising tendency from 2009 to 2018 (also from Period I to Period III ). Meng et al. [73] similarly concluded that during 2003 and 2014 in Beijing, SUHII exhibited a fluctuating increasing trend during summer days, and the same during the nights of all seasons, except for summer. However, the study of Liu et al. [43] revealed that the SUHII of Beijing during summer days showed a fluctuating increasing trend from 2003 to 2009, but from 2010 to 2018, it showed a fluctuating decreasing trend. The bias may be caused by different geographical focuses of Beijing and the calculation methods of SUHII.

The Combined Effects of ISA Change and the Greening Project on LST and SUHII
Through OLS estimation, it was found that afforestation and ISA changes together can only explain 11.59% to 22.47% of LST changes. However, previous studies have confirmed that the proportions of impervious surface (IS) and green spaces have a dominant effect on LST spatial differentiation [17,20]. The reason for this may be that the driving mechanism of LST changes in relation to LULC conversion is more complicated than with LULC patterns because of diversified conversion types and climate effects. OLS estimation may not be a perfect model for the studied variables, as it is based on a linear-based hypothesis and there is an assumption that the error terms are independent.
The results of OLS estimation showed that ΔISAKDE and ΔLST had a significant positive correlation on warm season days, and a significant negative correlation during cold season nights. In studies of pattern analysis, however, urban impervious surfaces and LST commonly show positive correlations during both the summer daytime/nighttime and winter nighttime [74,75]. Moreover, the results showed that FORKDE and ΔLST had a significant negative correlation on summer days, and the absolute value of the coefficient was much lower than that between ΔISAKDE and ΔLST. This showed a stronger role of ISA change of the same area in influencing LST on summer days, which is in accordance with previous studies [33][34][35]. During winter nights, the role of FORKDE in decreasing LST was similar to that of ΔISAKDE with the same area. This result matches the study of Duncan et al. [76], which showed that an increase in urban vegetation within a location reduces the LST in both summer and winter. However, it also challenges studies which concluded that increased greenery could reduce temperatures in summer, but increase them in winter [77,78]. The different results may be caused by the imparity of former land cover types at the revegetated sites, as well as the geographic and climatic context of studied areas.
The controlled regression analysis characterized the relationships between ΔLST, FORKDE, and ΔISAKDE more precisely, as it discussed the individual effects of FORKDE and ΔISAKDE on LST changes, respectively. When ΔISAKDE equaled zero, the relationship between ΔLST and FORKDE was similar to the result generated from OLS estimation, i.e., FORKDE decreased ΔLST both on warm season days and during cold season nights. A 10% increase of FORKDE resulted in a decrease of 0.53 °C of LST on summer days, and a decrease of 1.11 °C of LST during winter nights. This outcome is inconsistent with a previous study which concluded that vegetation cooling is generally stronger during the daytime periods, in warm seasons, in low-latitude zones, while in the opposite context, vegetation warming usually occurs [78]. The reason for this may be that most of the new plantations from the project were transformed from cropland; the effect of this conversion on influencing LST can be counteracted by the pre-existing role of previous land cover. When FORKDE equaled zero, the relationship between ΔLST and ΔISAKDE was different from the result of OLS estimation. On summer days and warm season days, a threshold of ΔISAKDE in relation to ΔLST was found, i.e., the magnitude of IS warming could be accelerated when ΔISAKDE exceeded 0.217 (about four times stronger than that when ΔISAKDE was below 0.217). During winter nights and cold season nights, a turning point of ΔISAKDE in relation to ΔLST was found (0.096 during winter nights and 0.034 during cold season nights), after which the negative correlation between ΔISAKDE and ΔLST reversed. This phenomenon implies that the effect of the IS dynamic on LST change is not monotonous. IS expansion heats the surface on warm season days, and the warming effect can be exacerbated after IS expansion reaches certain thresholds. During cold season nights, both IS shrinkage and strong IS expansion cause surface warming. The former could generally have been caused by land use conversion, but the latter was mainly caused by anthropogenic heat releases.
Seasonal variability of the LST changing pattern and SUHII was observed in this study, as well as in previous research [72,77,79,80], suggesting different driving mechanisms between year-round time divisions. When we assumed that LST variations during 2009-2018 in the plain area of Beijing were only caused by the combined effects of ISA changes and the planting project, it was found that 82.9% of the SUHII rise on warm season days (and 73.6% on summer days) was caused by the percent forest cover increase of 10.0, and 80.6% of the SUHII increase during cold season nights (and 78.9% during winter nights) could be attributed to the percent ISA increase of 6.6. This result conforms to the previous finding that SUHII on summer days is closely related to vegetation activity and anthropogenic heat releases, while during winter nights, SUHII is strongly linked to the albedo, builtup intensity (or night-time light), and anthropogenic heat releases [72,79]. However, for the first time, a quantitative relationship was demonstrated in this study. The mechanism behind SUHII changes may be as follows: On summer days and warm season days, large-scale afforestation in the periurban and rural areas lowered the regional ΔLST and enlarged the gap between the LST of urban and rural areas. IS expansion formed new urban areas and had a warming effect in these regions. Therefore, afforestation and ISA changes acted synergistically to increase SUHII. During winter nights and cold season nights, peri-urban and rural afforestation still decreased the regional ΔLST, thus increasing the gap between the LST of urban and rural areas. A small amount of slight ISA shrinkage in the high-urbanized areas increased the regional ΔLST, a large amount of medium ISA expansion in the suburban areas reduced the regional ΔLST, and a small amount of extensive ISA expansion in the rural areas (e.g., construction of the new Daxing Airport in the southern city periphery) increased the regional ΔLST. As a result, the scatter points of LST and ISAKDE in Period III showed a "U" shape trend compared with the predicted Scenario B, and the linear regression slope was eventually higher than that in the predicted scenario.

Implications for Future Landscape Planning for Urban Thermal Environment Regulation
This study characterized the combined effects of IS change and afforestation on LST and SUHII variations based on a spatial-temporal analysis of Beijing's plain area during 2009-2018. The results can provide important insights on UHI alleviation concerning IS and green space planning in a global context of rapid urbanization [29,67,81,82] and an anticipated new development stage when abundant urban (region) greening can be conducted to help cool the city [22,83]. First, large-scale afforestation in the periphery of built-up areas in Beijing was found to have a counterproductive effect on UHI attenuation. This was caused by a singular view of vegetation cooling, rather than considering the whole area and adapting the greening to the features of the IS distribution and its development. Therefore, it is recommended that IS and green spaces should be planned jointly in coping with the UHI effect.
Second, distinct diurnal and seasonal variations of LST and SUHII changes in response to ISA and forest cover dynamics were observed. On warm season days, LST change was more sensitive to the IS dynamic compared to vegetation change of the same area, while during winter nights, the difference was much smaller. SHUII changes in response to IS and vegetation cover variations seem to be more complex as the placing of IS and vegetation cover changes also matter. For example, in this study, the contribution of a 10.0% increase of vegetation cover versus a 6.6% increase of ISA to the SUHII increase was about 3:1 on summer days, but 1:4 during winter nights. Therefore, it is suggested that city planners are season-oriented when facing the task of UHI mitigation.
Third, the existence of threshold values of ISA change in influencing LST was identified. On warm season days, a greater regional ISA increase resulted in a higher LST increase, but when the regional ISA increase reached a certain threshold, the magnitude of the LST increase was exacerbated. During cold season nights, a greater regional ISA increase corresponded to a lower LST increase within a range, whereas after a threshold, a greater regional ISA increase resulted in a higher LST increase. As a result, practitioners and planners need to pay attention to the threshold of regional ISA expansion in relation to LST changes.

Limitations and Uncertainties
While our findings enrich the existing scientific knowledge of how IS and vegetation cover changes can influence the urban thermal environment, several study limitations should be highlighted. First, extra vegetation cover changes besides the afforestation project and the influence of other land cover conversions except for ISA and forest cover were omitted due to the restriction of appropriate data with a high accuracy. However, the error caused by this shortcoming is not believed to have had a major impact on the conclusions as no major forest gain/loss except for the planting project was recorded. The influence of additional land cover conversions was trivial, e.g., strict cropland protection was implemented by the government unless there was an official decision; bare soil and grassland only account for 9.7% and 1.1% of the study area, respectively; and the area of water cover was excluded during the calculation. Moreover, an in-depth analysis of urban thermal changes in response to various types of ISA and vegetated areas, and the mechanism of energy balance (natural or anthropogenic) behind this will be addressed in our future studies.
Second, inter-annual variations of meteorological elements (e.g., temperature, precipitation, and radiation level) in the same place may have influenced the results. However, those changes should be slight within a decade and given the city-wide scope. Additionally, the processes of the pixel group and data average during the calculation can reduce the spatial heterogeneity caused by this aspect. We found that the average surface LST in the whole city had a rising tendency during 2009-2018, probably partly caused by global climate change. Therefore, further studies can be conducted on the relationship between climate change and SUHII variation, but a longer time-series is needed.
Third, the revegetated areas in this study are still at an early growing stage, while the future effect of forest cover changes on LST may not be able to be interpreted in the same way. Therefore, results of longer observation periods and continuous research are needed. Moreover, the study only focused on the case of Beijing, whereas some of its conclusions may not be directly applicable elsewhere. Therefore, further work is needed in terms of global observations on the effect of LULC conversion (especially IS and vegetation cover variations) on LST and SUHII changes, in order to fully understand the characteristics of this relationship in different background climate and geographical contexts.

Conclusions
Knowledge of infrastructure-based interventions (especially ISA and green space regulations) to mitigate the UHI effect is imperative for planners and practitioners [7]. This study explored LST and SUHII changes in response to ISA and forest cover variations based on remote sensing analysis for the case of Beijing. The findings show that the percent of ISA and forest cover increased by 6.6 and 10.0, respectively, in the study area. SUHIIs showed significant rising tendencies during 2009-2018 in the four time divisions of warm season days, summer days, cold season nights, and winter nights. LST changes on warm season days and summer days responded positively to a regionalized ISA increase and negatively to a regionalized forest cover increase. However, during cold season nights and winter nights, LST changes responded negatively to a slight regional ISA increase, but positively to an extensive regional ISA increase, and LST variations responded negatively to a regional forest cover increase. The effect of vegetation cooling was weaker than ISA warming on warm season days, but the effect of vegetation cooling was similar to that of ISA during cold season nights. The afforestation project dominated the contribution of SUHII rise on warm season days, while ISA changes led the contribution of SUHII increase during cold season nights. Based on our results, we suggest that IS and green spaces should be planned jointly in coping with the UHI effect, strategies of UHI mitigation should be season-oriented, and the development threshold of ISA in influencing LST should be determined.