Effects of Green Space Patterns on Urban Thermal Environment at Multiple Spatial–Temporal Scales

Land use/land cover (LULC) pattern change due to human activity is one of the key components of regional and global climate change drivers. Urban green space plays a critical role in regulating urban thermal environment, and its cooling effect has received widespread attention in urban heat island (UHI) related studies. To fully understand the effects of the landscape pattern of an urban green space in regulating the urban thermal environment, it is necessary to further study the thermal effects of the landscape pattern of the urban green space and its characteristics under varied spatial–temporal scales. In this paper, we took the urban core area of Hangzhou City as the study area and analyzed the relationships between the landscape metrics of the urban green space and land surface temperature (LST) under varied spatial scales by using correlation analysis and redundancy analysis (RDA) methods. Multi-temporal Landsat 8 thermal infrared sensor data were used to retrieve the spatial and temporal dynamics of LSTs in four consecutive seasons, and the land use classification was interpreted using SPOT (Systeme Probatoire d’Observation de la Terre) satellite imagery. The results showed that landscape dominance metrics—e.g., percentage of landscape (PLAND) and largest patch index (LPI)—were the most influential factors on urban LST. The spatial configuration of urban landscape, as represented by the fragmentation and aggregation and connectedness, also showed significant effects on LST. Furthermore, landscape pattern metrics had varied spatial scale effects on LST. Specifically, the landscape dominance metrics of urban forest showed an increased influence on LST as the spatial scale increased, while for urban water, the trend was opposite. These findings might have some practical significance for urban planning about how to spatially arrange urban green space to alleviate UHI at local and regional scales.


Introduction
The land use/land cover (LULC) pattern change due to human activity is one of the key forces of regional and global climate change [1][2][3][4], and large-scale urbanization causing changes in land surface composition and structure during the past decades has been a main factor in LULC [5,6]. The expansion of urban land directly results in the relative reduction of natural vegetation and an increase of impervious surface. As a result, the changed evapotranspiration and thermal inertia [7][8][9] of the land surface unit, together with the heat produced by human activities during energy consumption [10,11],

Classification of Urban Land Use
Urban green spaces refer to all kinds of green infrastructures that use vegetated spaces and water bodies for the public well-being. Impervious surfaces are land surfaces, such as roads, parking lots, and building roofs, that repel rainwater and do not permit it to soak into the ground. It is noted that the study area is the core area of the city, so a large part of the land uses were built-up lands, mixed with urban forest composed by the mountainous area, wet land, parks, and trees along the roads. The cropland mainly distributed in the periphery of the urban core and its proportion was relatively small. We classified the land use in the study area into four Level I classes with 15 Level II sub-classes. The four Level I classes were Cropland (1), Forest (2), Water (3), and Built-up (4), and the corresponding 15 Level II sub-classes were Cropland (11), Street trees (21), Park trees (22), Mountain forest (23), Other forest (24), River (31), Lake (32), Pond (33), Residential (41), Industrial (42), Commercial (43), Administrative (44), Educational (45), Road (46), and Under construction (47). The first digit of the land classification code of the Level II class referred to its corresponding Level I class. The land use information in the study area was extracted from one cloudless SPOT 6 image obtained on 28 December, 2015. Imageries with higher spatial resolution could more accurately quantify the spatial pattern of urban land use, such as the smaller elements of greenspace, as well as changes in the configuration of larger patches [54]. The multispectral and panchromatic bands of the SPOT image were fused to get a 1.5 m resolution image using the Gram-Schmidt method. The SPOT image was radiometrically and atmospherically corrected, and the accuracy of the geometric correction was under 0.5 pixel. The pre-processed fused image was classified into 15 classes using the maximum likelihood classification method. The training sites of each Level II class were collected by field survey and visually identified polygons on the SPOT satellite imagery (R4, G3, and B2), and then each of the classes was assigned and recoded into four Level I classes, respectively ( Figure 2). The postclassification was assisted by using a 1:10,000 land use map and visual interpretation under the aid of Google Earth satellite images. The classification accuracy was assessed by stratified random sampling with reference to ground-truth data collected in the field survey. The assessment results showed that the overall accuracy of land use classification was above 95% and the kappa coefficient was 0.93.

Classification of Urban Land Use
Urban green spaces refer to all kinds of green infrastructures that use vegetated spaces and water bodies for the public well-being. Impervious surfaces are land surfaces, such as roads, parking lots, and building roofs, that repel rainwater and do not permit it to soak into the ground. It is noted that the study area is the core area of the city, so a large part of the land uses were built-up lands, mixed with urban forest composed by the mountainous area, wet land, parks, and trees along the roads. The cropland mainly distributed in the periphery of the urban core and its proportion was relatively small. We classified the land use in the study area into four Level I classes with 15 Level II sub-classes. The four Level I classes were Cropland (1), Forest (2), Water (3), and Built-up (4), and the corresponding 15 Level II sub-classes were Cropland (11), Street trees (21), Park trees (22), Mountain forest (23), Other forest (24), River (31), Lake (32), Pond (33), Residential (41), Industrial (42), Commercial (43), Administrative (44), Educational (45), Road (46), and Under construction (47). The first digit of the land classification code of the Level II class referred to its corresponding Level I class. The land use information in the study area was extracted from one cloudless SPOT 6 image obtained on 28 December, 2015. Imageries with higher spatial resolution could more accurately quantify the spatial pattern of urban land use, such as the smaller elements of greenspace, as well as changes in the configuration of larger patches [54]. The multispectral and panchromatic bands of the SPOT image were fused to get a 1.5 m resolution image using the Gram-Schmidt method. The SPOT image was radiometrically and atmospherically corrected, and the accuracy of the geometric correction was under 0.5 pixel. The pre-processed fused image was classified into 15 classes using the maximum likelihood classification method. The training sites of each Level II class were collected by field survey and visually identified polygons on the SPOT satellite imagery (R4, G3, and B2), and then each of the classes was assigned and recoded into four Level I classes, respectively ( Figure 2). The post-classification was assisted by using a 1:10,000 land use map and visual interpretation under the aid of Google Earth satellite images. The classification accuracy was assessed by stratified random sampling with reference Sustainability 2020, 12, 6850 4 of 18 to ground-truth data collected in the field survey. The assessment results showed that the overall accuracy of land use classification was above 95% and the kappa coefficient was 0.93.

LST Retrieval
The Landsat satellite series with a thermal infrared band is the most common imagery source for LST retrieval. From the Landsat 4 Thematic Mapper to Landsat 8, several widely used algorithms have been developed for retrieving LST from the thermal channel of Landsat imageries: for example, the radiative transfer equation (RTE) and single channel algorithm (SCA) [55,56], the mono-window algorithm (MWA) [57,58], and the spit-window algorithm (SWA) [59].
LST observations acquired by remote sensing technologies have been used to assess the UHI, to develop land surface-atmosphere exchange models, and to analyze the relationship between temperature and land use and land cover in urban areas. The multi-temporal satellite data is effective and useful in diurnal, seasonal, and interannual urban land surface temperature variation studies [12,60]. In this study, we retrieved LSTs from four Landsat 8 imageries representing four consecutive seasons. The acquisition dates of the Landsat 8 images were 13 October, 2015, 16 December, 2015, 15 March, 2016, and 27 July, 2016, all at 10:31 a.m. local time. The 13 October and 16 December represented autumn and winter, the cold and dry seasons, respectively; the 15 March represented spring, the dry and warm season; and 27 July represented summer, the hot and humid season. The images were generally cloud free and the imaging dates were approximately in the middle of the seasons. These dates were no more than 8 months from the acquisition date of the SPOT image used for urban green space mapping, and it was reasonable to assume that there was no significant change in urban land uses [45]. To retrieve LST, Landsat 8 Thermal Infrared Sensor (TIRS) band 10 is more commonly used than band 11 due to its higher accuracy [61,62], and the RTE applied is: where Lsensor is the radiance measured by the sensor, ε is the land surface emissivity, Ts is the land surface temperature according to B(Ts), which is the blackbody radiance given by the Planck's law, τ is the total atmospheric transmissivity for channel λ between the surface and the sensor, and Lup and Ldown, are the downwelling and upwelling atmospheric radiance, respectively. The atmospheric parameters τ, Lup and Ldown can be obtained by MODTRAN (MODerate resolution atmospheric TRANsmission) via the atmospheric correction parameter calculator provided by NASA (https://atmcorr.gsfc.nasa.gov/) [63]. The value of Lsensor can be directly converted from Landsat 8 TIRS band data by using the radiance rescaling factors provided in the metadata file

LST Retrieval
The Landsat satellite series with a thermal infrared band is the most common imagery source for LST retrieval. From the Landsat 4 Thematic Mapper to Landsat 8, several widely used algorithms have been developed for retrieving LST from the thermal channel of Landsat imageries: for example, the radiative transfer equation (RTE) and single channel algorithm (SCA) [55,56], the mono-window algorithm (MWA) [57,58], and the spit-window algorithm (SWA) [59].
LST observations acquired by remote sensing technologies have been used to assess the UHI, to develop land surface-atmosphere exchange models, and to analyze the relationship between temperature and land use and land cover in urban areas. The multi-temporal satellite data is effective and useful in diurnal, seasonal, and interannual urban land surface temperature variation studies [12,60]. In this study, we retrieved LSTs from four Landsat 8 imageries representing four consecutive seasons. The acquisition dates of the Landsat 8 images were 13 October, 2015, 16 December, 2015, 15 March, 2016, and 27 July, 2016, all at 10:31 a.m. local time. The 13 October and 16 December represented autumn and winter, the cold and dry seasons, respectively; the 15 March represented spring, the dry and warm season; and 27 July represented summer, the hot and humid season. The images were generally cloud free and the imaging dates were approximately in the middle of the seasons. These dates were no more than 8 months from the acquisition date of the SPOT image used for urban green space mapping, and it was reasonable to assume that there was no significant change in urban land uses [45]. To retrieve LST, Landsat 8 Thermal Infrared Sensor (TIRS) band 10 is more commonly used than band 11 due to its higher accuracy [61,62], and the RTE applied is: where L sensor is the radiance measured by the sensor, ε is the land surface emissivity, T s is the land surface temperature according to B(T s ), which is the blackbody radiance given by the Planck's law, τ is the total atmospheric transmissivity for channel λ between the surface and the sensor, and L up and L down , are the downwelling and upwelling atmospheric radiance, respectively. The atmospheric parameters τ, L up and L down can be obtained by MODTRAN (MODerate resolution atmospheric TRANsmission) via the atmospheric correction parameter calculator provided by NASA (https://atmcorr.gsfc.nasa.gov/) [63]. The value of L sensor can be directly converted from Landsat 8 TIRS band data by using the radiance rescaling factors provided in the metadata file of Landsat 8 data [62]. The value of ε was estimated by the NDVI thresholds method (NTM), which uses certain NDVI values to distinguish between soil pixels, full vegetation pixels, and mixed pixels composed of soil and vegetation [64]. The emissivities of water and built-up areas were adopted from the ASTER Global Emissivity Database 4.0 [58,65].
With the thermal radiance measured at sensor level, the atmospheric parameters obtained from the radiative transfer model and the land surface emissivity estimated using NTM, the T s could be retrieved as follows [61]: where K 1 is 774.89 Wm −2 sr −1 µm −1 , and K 2 is 1321.08 K [66].

Landscape Pattern Metrics
The landscape metrics are quantitative representations of the spatial composition of a particular landscape pattern. In this study, we used four categories of landscape pattern metrics-i.e., dominance, shape complexity, aggregation and connectedness, and fragmentation [39][40][41][42][43][44]-to quantify the landscape pattern and explore the corresponding correlations with LST (Table 1 and Table S1 in the supplementary materials). The landscape pattern metrics were analyzed by the software FRAGSTATS 4.3.

Statistical Analysis
To explore the effects of urban green space patterns on LST, different grid scales were utilized. In the literature, the grid size generally ranges from 200 to 1000 m, or even 3000 m, according to the greenspace condition or the sampling density in a case-specific manner [40,42,68,69]. In this study, we repeatedly divided the study area into fishnet-like grids with different spatial scales ( Figure 3)-i.e., 0.5 to 4.0 km every 500 m [44]-to study the relationships between the landscape metrics of the urban green space and LST under varied spatial scales. Since this study mainly concerned the influence of the landscape pattern of urban green space on LST, the sampling grids containing cropland parcels were excluded; in addition, the grids at the edge of the study area with blank areas of more than 20% of the total area of the specific grid were also excluded, including the grids with cloud or cloud shadow. As a result, the valid sampling grid numbers for the spatial scale ranging from 0.5-4.0 km every 500 m were 2000, 394, 156, 77, 41, 29, 23, and 19, respectively. Sustainability 2020, 12, x FOR PEER REVIEW 6 of 18

Statistical Analysis
To explore the effects of urban green space patterns on LST, different grid scales were utilized. In the literature, the grid size generally ranges from 200 to 1000 m, or even 3000 m, according to the greenspace condition or the sampling density in a case-specific manner [40,42,68,69]. In this study, we repeatedly divided the study area into fishnet-like grids with different spatial scales ( Figure 3)i.e., 0.5 to 4.0 km every 500 m [44]-to study the relationships between the landscape metrics of the urban green space and LST under varied spatial scales. Since this study mainly concerned the influence of the landscape pattern of urban green space on LST, the sampling grids containing cropland parcels were excluded; in addition, the grids at the edge of the study area with blank areas of more than 20% of the total area of the specific grid were also excluded, including the grids with cloud or cloud shadow. As a result, the valid sampling grid numbers for the spatial scale ranging from 0.5-4.0 km every 500 m were 2000, 394, 156, 77, 41, 29, 23, and 19, respectively. The Pearson correlation analysis and direct gradient analysis were used to analyze the relationships between urban LST and the landscape metrics of the urban green space. In environmental studies, one of the most common tasks is to identify the most influential environmental factors. Compared with a regression analysis, factor analysis, or principal analysis, the direct gradient analysis is special in its ability to quantify the unique variation belonging to specific explanatory variables [70][71][72]. In the direct gradient analysis, the detrended correspondence analysis (DCA) was used to detect gradient length in the data sets. In the case of a short gradienti.e., the longest gradient length obtained with DCA was less than 3 standard deviations-the linear model with redundancy analysis (RDA) was applied; otherwise, the canonical correspondence analysis (CCA) was used. We applied DCA to LST and found that all gradient lengths were less than 3 standard deviations, so RDA was applied to explore the relationships among LST and landscape metrics: i.e., the species and environmental variables in the RDA context, respectively. Two important outputs could be obtained by RDA: i.e., the proportions (%) of the total variance of LST explained by each canonical axis, and the ordination diagrams, which explicitly reveal the relationships between landscape pattern metrics and LST. A positive correlation is expected when the arrows of two variables point to the same direction, and vice versa. The RDAs were performed using the CANOCO 4.5 program [73]. The Pearson correlation analysis and direct gradient analysis were used to analyze the relationships between urban LST and the landscape metrics of the urban green space. In environmental studies, one of the most common tasks is to identify the most influential environmental factors. Compared with a regression analysis, factor analysis, or principal analysis, the direct gradient analysis is special in its ability to quantify the unique variation belonging to specific explanatory variables [70][71][72]. In the direct gradient analysis, the detrended correspondence analysis (DCA) was used to detect gradient length in the data sets. In the case of a short gradient-i.e., the longest gradient length obtained with DCA was less than 3 standard deviations-the linear model with redundancy analysis (RDA) was applied; otherwise, the canonical correspondence analysis (CCA) was used. We applied DCA to LST and found that all gradient lengths were less than 3 standard deviations, so RDA was applied to explore the relationships among LST and landscape metrics: i.e., the species and environmental variables in the RDA context, respectively. Two important outputs could be obtained by RDA: i.e., the proportions (%) of the total variance of LST explained by each canonical axis, and the ordination diagrams, which explicitly reveal the relationships between landscape pattern metrics and LST. A positive correlation is expected when the arrows of two variables point to the same direction, and vice versa. The RDAs were performed using the CANOCO 4.5 program [73].
A comprehensive flowchart was designed to show the variables used in the statistical analysis and how they were combined with the spatial and temporal analysis (Figure 4). Firstly, the satellite imagery-based land use and LST data were retrieved from SPOT 6 and Landsat 8 images, respectively, and then the landscape pattern metrics and mean LST of each grid at different spatial scales were calculated under varied grid scales ranging from 0.5 to 4 km. Finally, the relationships between urban green space patterns and LST at multiple spatial-temporal scales were analyzed by using a correlation analysis and RDA. A comprehensive flowchart was designed to show the variables used in the statistical analysis and how they were combined with the spatial and temporal analysis (Figure 4). Firstly, the satellite imagery-based land use and LST data were retrieved from SPOT 6 and Landsat 8 images, respectively, and then the landscape pattern metrics and mean LST of each grid at different spatial scales were calculated under varied grid scales ranging from 0.5 to 4 km. Finally, the relationships between urban green space patterns and LST at multiple spatial-temporal scales were analyzed by using a correlation analysis and RDA.

Spatial Pattern of LST across Different Seasons
There were remarkable seasonal and spatial variations in LST. Figure 5 shows the spatial distribution of LSTs in four seasons. The average LSTs were 24.93, 10.42, 25.71, and 34.58 °C in autumn (October), winter (December), spring (March), and summer (July), respectively ( Table 2). The maximum standard deviation (SD) of LST, 4.23 °C, was in July, followed by 3.33, 2.37, and 1.32 °C in March, October, and December, respectively.

Spatial Pattern of LST across Different Seasons
There were remarkable seasonal and spatial variations in LST. Figure 5 shows the spatial distribution of LSTs in four seasons. The average LSTs were 24.93, 10.42, 25.71, and 34.58 • C in autumn (October), winter (December), spring (March), and summer (July), respectively ( Table 2). The maximum standard deviation (SD) of LST, 4.23 • C, was in July, followed by 3.33, 2.37, and 1.32 • C in March, October, and December, respectively.  Table 3 shows the LST statistics of the different land uses in four seasons. The impervious surface generally had the highest mean LSTs in four seasons, except in December, 2015, when its mean LST was equal to the LST of urban water. The urban water had the lowest mean LSTs in all seasons except winter (December). The mean LSTs of urban forest were higher than those of urban water (except December), but lower than cropland. The largest LST variation was found in summer (July) when the impervious surface and urban forest had the largest temperature difference compared with the other seasons. In contrast, the mean LSTs of all the land uses were much closer in winter (December).   Table 3 shows the LST statistics of the different land uses in four seasons. The impervious surface generally had the highest mean LSTs in four seasons, except in December, 2015, when its mean LST was equal to the LST of urban water. The urban water had the lowest mean LSTs in all seasons except winter (December). The mean LSTs of urban forest were higher than those of urban water (except December), but lower than cropland. The largest LST variation was found in summer (July) when the impervious surface and urban forest had the largest temperature difference compared with the other seasons. In contrast, the mean LSTs of all the land uses were much closer in winter (December).

Relationships between Land Use Proportion and LST
The correlation analysis showed that the proportions of vegetation and water at grid level had significant negative correlations with LST on most spatiotemporal scales, except winter for water; for the impervious surface, significant positive correlations with LST could be found on all scales in each season ( Figure 6). The absolute value of the correlation coefficients between forest proportion and LST was the largest in summer, followed by autumn, spring, and winter, and increased with the increase of grid scales. The water proportion was significantly negatively correlated with LST in spring, summer, and autumn, but significantly positively correlated with LST in winter at small scales: e.g., 0.5-2.0 km. The absolute values of the correlation coefficients decreased with the increase of grid scales. The positive correlation coefficients between impervious surface and LST increased with the increase of grid scales, and were larger in summer and autumn, followed by spring and winter.
increase of grid scales. The water proportion was significantly negatively correlated with LST in spring, summer, and autumn, but significantly positively correlated with LST in winter at small scales: e.g., 0.5-2.0 km. The absolute values of the correlation coefficients decreased with the increase of grid scales. The positive correlation coefficients between impervious surface and LST increased with the increase of grid scales, and were larger in summer and autumn, followed by spring and winter.

Relationships between Urban Green Space Landscape Pattern and LST
In the Pearson correlation analysis, the correlation coefficients between land uses proportion (at grid scale) and LST increased with the increasing of the spatial scale ( Figure 6). To simplify the analysis of the effects of landscape patterns on the spatial variation of urban LSTs, RDA was applied

Relationships between Urban Green Space Landscape Pattern and LST
In the Pearson correlation analysis, the correlation coefficients between land uses proportion (at grid scale) and LST increased with the increasing of the spatial scale ( Figure 6). To simplify the analysis of the effects of landscape patterns on the spatial variation of urban LSTs, RDA was applied at 1, 2, and 4 km grid scales, respectively. A Monte Carlo permutation test (499 permutations) was used to determine the statistical validity of RDA, and the significant explanatory variables were identified. The explanations of landscape pattern metrics on LST at different spatial scales and the biplots of LST and significant landscape metrics according to RDA were shown in Table 4 and Figure 7, respectively.

Seasonal Variability in the Relationships between Urban Green Space and LST
In this study, the mean LSTs of impervious surface in all seasons except winter (December) were higher than the other land uses (Table 3), and the correlation analysis further indicated that LST was positively correlated with the proportion of impervious surface in all seasons (Figure 6), while negatively correlated with that of urban forest; for urban water, the correlation with LST was negative except in winter. This was consistent with the general consensus that vegetation and water could alleviate UHI, while impervious surface was the opposite due to the different specific heat capacities [74][75][76][77][78]. Furthermore, we showed that the mean LSTs of different land uses differed mostly in July, when the mean LSTs of the impervious surface were 4.4 and 8.7 °C higher than the urban forest and urban water, respectively. The mean LSTs of different land uses were closest in December: for Figure 7. Biplots of LST (represented by blue lines) and landscape metrics (red lines) according to RDA. RDA was conducted separately under 1 km (a-c), 2 km (d-f), and 4 km (g-i) scales. The length of the arrow represents the standard deviation of variables in the sorting space, and the direction indicates the change of gradient direction. The cosine value of the angle between the arrow and axis represents the correlation. The explanatory landscape pattern variables were selected according to p-value < 0.05. At 1, 2, and 4 km grid scales, the 32.9%, 63.7%, and 99.8% variations in LST could be explained by the landscape patterns of urban forest, respectively; the corresponding percentages for urban water were 71.6%, 82.8%, and 99.3%, respectively; for the impervious surface, the percentages were 69.7%, 84.6%, and 98.8%, respectively. The landscape metrics of PLAND, LPI, area-weighted mean perimeter-area ratio (PAAM), landscape division index (DIVISION), PLAND, AI, cohesion index (COHESION), and area-weighted mean contiguity index (CONAM) showed relatively high explanatory capabilities at all three scales, in which PLAND, LPI, and DIVISION were the most influential explanatory variables. There were 4%, 8%, and 11% variations in LST at 1, 2, and 4 km grid scales, which could be explained by the landscape dominance metrics PLAND and LPI, respectively; for urban water, about 8% of variations in LST could be explained by PLAND and LPI at all three scales, but the explained variations slightly decreased as the spatial scale increased; for the impervious surface, all the explained variations in LST by PLAND and LPI were greater than 10% but decreased as the spatial scale increased. The metric DIVISION, which represented the landscape fragmentation of urban forest, urban water, and impervious surface, showed similar explanatory capabilities to LST variations, such as PLAND and LPI, at all three scales.
The landscape complexity metrics-except mean perimeter-area ratio (PAMN) and LSI-of urban water showed relative high explanatory capabilities on LST compared with those of urban forest and impervious surface, and the explained variations in LST increased as the spatial scale increased (Table 4). For example, the explained variations in LST by mean shape index (SHMN), area-weighted mean shape index (SHAM), mean fractal dimension index (FRMN), and area-weighted mean fractal dimension index (FRAM) of the urban water increased from 3.6%-6.7% at 1 km to 6.1%-8.1% and 7.3%-13.4% at 2 and 4 km spatial scales, respectively. For the aggregation and connectedness metrics, the explanatory capabilities of proportion of like adjacencies (PLADJ), AI, COHESION, and CONAM of the urban forest and impervious surface increased as the spatial scale increased; the explained variations in LST by COHESION of urban water showed an increasing trend as the spatial scale increased, while PLADJ, AI, and CONAM were the opposite. Figure 7 shows the ordination diagrams of the RDA analysis of seasonal LST with respect to the landscape metrics of urban forest, urban water, and impervious surface at 1, 2, and 4 km grid scales, respectively. The first axis in the ordination diagrams represents changes of LST in three seasons (March, July, and October), and in one season (December) to a lesser extent. The arrows of two variables pointing to the same direction indicate a positive correlation, and the angle between two arrows is inversely proportional to the degree of their correlation. The LSTs in all seasons were negatively correlated with the landscape dominance, and the aggregation and connectedness metrics-e.g., the PLAND, LPI, PLADJ, AI, and clumpiness index (CLUMPY)-of urban forest at all scales, and were positively correlated with landscape fragmentation and complexity metrics: e.g., DIVISION and PAAM, at all scales. The correlations between LST and the landscape metrics of urban water were not consistent among seasons. In March, July, and October, urban water had similar effects on LST as those of urban forest; however, the correlations between LST in December and landscape dominance, aggregation and connectedness, fragmentation, and shape complexity-as expressed by the metrics of PLAND, LPI, PLADJ, AI, CLUMPY, DIVISION, and PAAM-of urban water were much smaller compared with the other seasons. The impervious surface showed opposite effects on LST compared with that of urban forest. The LSTs in all seasons were negatively correlated with the landscape shape complexity and fragmentation indices-e.g., PAAM and DIVISION-at all spatial scales, and negatively correlated with landscape dominance and aggregation and connectedness of the impervious surface: e.g., PLAND, LPI, and PLADJ.

Seasonal Variability in the Relationships between Urban Green Space and LST
In this study, the mean LSTs of impervious surface in all seasons except winter (December) were higher than the other land uses (Table 3), and the correlation analysis further indicated that LST was positively correlated with the proportion of impervious surface in all seasons (Figure 6), while negatively correlated with that of urban forest; for urban water, the correlation with LST was negative except in winter. This was consistent with the general consensus that vegetation and water could alleviate UHI, while impervious surface was the opposite due to the different specific heat capacities [74][75][76][77][78]. Furthermore, we showed that the mean LSTs of different land uses differed mostly in July, when the mean LSTs of the impervious surface were 4.4 and 8.7 • C higher than the urban forest and urban water, respectively. The mean LSTs of different land uses were closest in December: for example, the LST difference between urban forest and impervious surface was only 1.3 • C; the LST differences between impervious surface and the other land uses in October and March were in between those of July and December (Table 3). This indicated that different land uses might have varied thermal environment effects.
Numerous studies have showed that there are seasonal differences in the cooling effect of urban green space-generally, the strongest is in summer while the weakest is in winter-and it might have a close relationship with the evapotranspiration of the urban forest [79,80]. This was consistent with the Pearson correlation analysis between LSTs and the proportion of urban forest on multiple spatial scales ( Figure 6), in which the correlation coefficients in warmer seasons were much higher than those in winter. Previous studies also showed that the correlations between land uses and LST in winter were relatively weak compared with the other seasons [81,82]. The RAD analysis also showed that landscape pattern metrics could explain more variations in LST in July than December (Figure 7). This suggested that integrating urban forests with urban water, such as urban rivers, lakes, and ponds, may have great potential to improve human thermal comfort in urban areas. Additionally, although cropland is slightly weaker than that of urban water in regulating the thermal environment all year around, it should be strengthened in the future together with the optimal allocation and design of other urban green spaces as a sustainable food production system with economic and environmental benefits.

Influential Landscape Pattern Metrics on LST
The RDA analysis showed that landscape dominance metrics PLAND and LPI were the most influential factors on urban LST (Table 4). Generally, the explanatory capabilities of the landscape dominance of urban forest and water were higher than the other landscape pattern metrics, and this suggested that the cooling effects might be closely related with the enhanced evapotranspiration as a result of increased forest and water coverage at local scales. In the literature, the land use composition was found to be more influential on LST than its spatial configuration [41,43,83,84], and this was consistent with our findings. This means that increasing vegetation coverage and expanding the water area in urban areas are still the effective measures to increase surface cooling coverage.
However, we found that the spatial configuration of the urban landscape also showed significant effects on LST. Specifically, the explanatory capability of DIVISION, one metric of the landscape fragmentation, was almost comparable with that of LPI, and it meant that the more fragmented and dispersed the urban forest, including water area and contiguous sprawl of construction land, the higher the LST. Meanwhile, the PD and number of patches (NP) of urban forest showed significant positive correlations with LST, and this implied that the fragmentation of urban green space might be an important factor leading to the increase of LST; similar reports could be found in [41][42][43]85]. On the contrary, compared with DIVISION, the more the aggregation and connectedness of the urban forest and urban water, the better the cooling effect, while for the impervious surface, the trend was opposite.
In terms of implications for urban planning and land use management, it might still be necessary to increase the proportion of urban green spaces-just as with the common practices in urban green space planning-to relieve UHI. However, the focus should also be on optimizing the layout of urban green spaces to reduce LST due to the shortage of urban land [30,68,86]. The relationships of landscape fragmentation and aggregation and connectedness of urban green space to LST revealed in this study means that, under the premise that the area of urban green space remains unchanged, the larger and more concentrated and connected urban green spaces might be more effective in cooling land surface, rather than numerous small green space patches.

Scaling Effects of Urban Green Space on LST
The landscape pattern metrics could explain more than 60% variation in the urban thermal environment at all spatial scales except at the 1 km grid scale, in which 32.9% of variations in LST were explained (Table 4). This confirmed that the composition and configuration of urban landscape might be the key factors controlling the spatial variation of LST [68,69].
Meanwhile, the landscape pattern metrics showed varied spatial scale effects on LST. For example, the RDA analysis showed that the landscape dominance metrics (PLAND and LPI) of urban forest had an increased influence on LST as the spatial scale increased, while for urban water, the trend was opposite (Figure 7). The correlation analysis also revealed that the proportions of urban forest and urban water had opposite trends in correlation with LST as the spatial scale increased ( Figure 6). Additionally, we found that a large part of landscape metrics-for example, the aggregation and connectedness ones of urban forest, urban water, and impervious surface, and the shape complexity ones of urban water-had increased effects on LST as the spatial scale increased. These findings might have some practical significance for urban planning about how to arrange urban green space to alleviate UHI. On a smaller block scale, for example, a reasonably increasing urban water area, if possible, and avoiding the continuous development of construction land might achieve a certain cooling effect, while on urban or regional scales, it might be more effective to alleviate UHI by optimizing the spatial pattern of urban forest, such as increasing the coverage and improving the agglomeration and connectivity of urban forest. Planning and managing urban green space and urban heat island mitigation needs to be approached holistically, taking into account diverse spatial-temporal dynamics [87,88].

Conclusions
There were remarkable seasonal and spatial variations in LST correlated with specific land uses. The impervious surface and urban water had the highest and lowest mean LSTs in all seasons except winter (December), respectively. The largest LST variation was found in summer (July) compared with the other seasons. In contrast, the mean LSTs of all the land uses were much closer in winter (December). In all seasons, the LSTs were positively correlated with the proportion of impervious surface, while negatively correlated with that of urban forest; for urban water, the correlation was negative, except in winter. The cooling effect of urban forest showed seasonal differences: generally, the strongest were in summer while the weakest were in winter.
Generally, the explanatory capabilities of landscape dominance metrics were higher than the other landscape pattern metrics. The spatial configuration of urban landscape, as represented by the fragmentation and aggregation and connectedness of the land use patches, also showed significant effects on LST. We suggested herein that it might still be necessary to increase the proportion of urban green spaces, just as with the common practices applied in urban green space planning, to relieve UHI; however, for cases where there is no more land allocated for urban green space, an optimized spatial arrangement of urban green space might also be an available means to counteract UHI.
The landscape pattern metrics had varied spatial scale effects on LST. Specifically, the landscape dominance metrics of urban forest showed that the influence on LST increased as the spatial scale increased, while for urban water, the trend was opposite. These findings might have some practical significance for urban planning about how to spatially arrange urban green space to alleviate UHI at local and regional scales.
Although apparently conclusive, there were limitations in interpreting the relationship between landscape pattern metrics and LST at multiple spatial scales. Numerous studies about the relationship between LST and urban green space are still case-specific. Further research should also consider the influential factors of urban morphology, such as the tree canopy cover, urban block form, and mountain topography, and comparison studies across different cities under varied climatic conditions will be of particular interest. Additionally, it would be helpful to establish a common set of landscape metrics to avoid redundancy and improve the comparability across similar studies.
Author Contributions: Y.S. and X.S. conceived and designed the framework of the research structure. All authors contributed to the data analysis, interpretation, and manuscript writing. All authors have read and agreed to the published version of the manuscript.