The Effect of Urban Green Spaces on the Urban Thermal Environment and Its Seasonal Variations

Urban green spaces have been shown to decrease land surface temperature (LST) significantly. However, few studies have explored the relationships between urban green spaces and LST across different seasons at different spatial scales. In this study, using Changchun, China as a case study, landscape ecology and comparative approaches were employed quantitatively to investigate the effects of the composition and configuration of urban green spaces on the urban thermal environments. LST maps were retrieved from Landsat 8 Thermal Infrared Sensor (TIRS) data acquired on four dates that represented four different seasons, and detailed information of urban green spaces was extracted from high resolution imagery GF-1. Normalized differential vegetation index (NDVI) and six landscape metrics at patch, class, and landscape level were used to characterize the spatial patterns of urban green spaces. The results showed that urban green spaces did have significant cooling effects in all seasons, except for winter, but the effects varied considerably across the different seasons and green types, and seemed to depend on the NDVI and size of urban green spaces. Compared to shape metrics, the negative relationships between the LST and the area and the NDVI of urban green spaces were more significant. Both the composition and configuration of urban green spaces can affect the distribution of LST. Based on findings with one city, given a fixed area of urban green spaces, the number of green patches can positively or negatively affect the LST, depending on if the number is larger than a threshold or not, and the threshold varies according to the given area. These findings provide new perspectives, and further research is also suggested, to generate a better understanding of how urban green spaces affect the urban thermal environment.


Introduction
According to the United Nations, over half of the world's population (54%) lived in urban areas in 2014, and this number is projected to increase to 66% by 2050 as a result of rapid urbanization [1]. Therefore, more people will suffer from urban heat island (UHI) effects, which are caused by dramatic changes in land use and land cover, energy use, and air-pollutant emissions [2][3][4]. Intensive heat can cause many negative effects, including significant health and well-being problems for urban inhabitants over the long term, and especially for the elderly [5][6][7]. There is no doubt that understanding how spaces at multiple different spatial scales rather than just one single scale which is the spatial resolution of the imagery data.
This paper aims to identify appropriate metrics that effectively describe the features of urban green spaces in order to investigate how different green spaces and their spatial configurations affect the urban thermal environments at different spatial scales across different seasons and to establish quantitative relationships between these aspects. To achieve our goals, very high spatial resolution imagery GF-1 (2-m spatial resolution) was utilized to obtain detailed information of urban green spaces, and Landsat 8 TIRS (30-m spatial resolution) was chosen to provide data on the urban thermal environment in spring, summer, autumn, and winter in Changchun, China.

Study Area
Changchun city (125 • 06 -125 • 36 E, 43 • 43 -44 • 04 N) (Figure 1), which is the capital of Jilin province, China, is located in the northern hemisphere mid-latitude regions, with an urban population of 3.65 million (2014). The city serves as the financial center and cultural heart of the province. In the past 30 years, Changchun has been undergoing an accelerated process of urbanization, with the urban built-up area having increased from 138.4 km 2 in 1984 to 557.9 km 2 in 2014 [44]. According to the National Economy and Society statistical bulletin of Changchun in 2014, the proportion of urban green area in the Changchun built-up area is approximately 36.5%, and the per capita green area is about 11.6 m 2 . This paper aims to identify appropriate metrics that effectively describe the features of urban green spaces in order to investigate how different green spaces and their spatial configurations affect the urban thermal environments at different spatial scales across different seasons and to establish quantitative relationships between these aspects. To achieve our goals, very high spatial resolution imagery GF-1 (2-m spatial resolution) was utilized to obtain detailed information of urban green spaces, and Landsat 8 TIRS (30-m spatial resolution) was chosen to provide data on the urban thermal environment in spring, summer, autumn, and winter in Changchun, China.

Study Area
Changchun city (125°06′-125°36′ E, 43°43′-44°04′ N) (Figure 1), which is the capital of Jilin province, China, is located in the northern hemisphere mid-latitude regions, with an urban population of 3.65 million (2014). The city serves as the financial center and cultural heart of the province. In the past 30 years, Changchun has been undergoing an accelerated process of urbanization, with the urban built-up area having increased from 138.4 km 2 in 1984 to 557.9 km 2 in 2014 [44]. According to the National Economy and Society statistical bulletin of Changchun in 2014, the proportion of urban green area in the Changchun built-up area is approximately 36.5%, and the per capita green area is about 11.6 m 2 . The region is characterized by a sub-humid continental climate, four seasons with a long, cold winter and a short, hot summer. According to meteorological records of 1951-2014, the annual precipitation of Changchun is 561.6 mm and the annual average temperature is 5.5 °C, with a hot July (23.1 °C) and a cold December (−12.3 °C). The annual average air temperature has increased by 1.86 °C during last 60 years from 1951 to 2011 [38]. As a result of urbanization, the average urban air temperature is about 0.1-0.5 °C higher than that of surrounding rural areas based on the local meteorological data.
There are few studies using Landsat 8 data and very high resolution images to estimate the relationships between urban green spaces and thermal environments in the cold temperate zone. This study may help to provide new insights for UHI mitigation and urban planning. The region is characterized by a sub-humid continental climate, four seasons with a long, cold winter and a short, hot summer. According to meteorological records of 1951-2014, the annual precipitation of Changchun is 561.6 mm and the annual average temperature is 5.5 • C, with a hot July (23.1 • C) and a cold December (−12.3 • C). The annual average air temperature has increased by 1.86 • C during last 60 years from 1951 to 2011 [38]. As a result of urbanization, the average urban air temperature is about 0.1-0.5 • C higher than that of surrounding rural areas based on the local meteorological data. There are few studies using Landsat 8 data and very high resolution images to estimate the relationships between urban green spaces and thermal environments in the cold temperate zone. This study may help to provide new insights for UHI mitigation and urban planning.

Urban Green Spaces Extraction
Detailed information of urban green spaces was derived from a GF-1 image acquired on 22 June 2015 that was provided by the China Center for Resources Satellite Data and Application (CRESDA). GF-1 is equipped with two multispectral scanners with 2 m resolution panchromatic and 8 m resolution multispectral, and 4 multispectral scanners with 16 m resolution. Satellite engineering made a breakthrough in optical remote sensing technology by integrating high space resolution, multispectrum, and high time resolution. The multispectural and panchromatic images were first fused to produce a four-band pan multispectral image with a resolution of 2 m. The images fusion was performed in ENVI 5.3. First, RPC Orthorectification Workflow was applied to the multispectral (4 m) and panchromatic (2 m) images, then Image Registration Workflow was used to perform image registration. Finally, the NNDiffuse Pan Sharpening method was used to fuse images. The urban green spaces were manually classified into six types based on their locations and services, as defined in Table 1. Visual interpretation was used to extract information about the green spaces, including the boundary and type of each patch, according to location, size, shape, and spatial relationships with neighbors. In addition, Google Earth maps with higher spatial resolution were combined to identify the type of patch. In our case, the minimum patch size is 200 m 2 . Fieldwork conducted during August 2015 resulted in the acquisition of 96 ground references for accuracy assessment. Overall, the total accuracy of the derived urban green map achieved 92.1%. Table 1. Urban green spaces types and its descriptions.

Types and Abbreviation Descriptions
Attached Green Spaces, AGS green spaces located in commercial , residential areas and institutions Road Green Spaces, RGS green spaces on the sides of road, usually shown as belt bands Park Green Spaces, PGS theme parks, comprehensive parks, forest parks, and other parks Productive Plantation Green Spaces, PPGS nursery gardens, featuring saplings, flowers, and plants Ecological Green Spaces, EGS windbreaks, sediment control, used for environment protection Other Green Spaces, OGS other green spaces types not belonging to any of the above types

Land Surface Temperature
The successful launch of Landsat 8 on 11 February 2013 with the most advanced Thermal Infrared Sensor (TIRS) is anticipated to extend the 40-year Landsat record for at least another 5 years [31]. A series of Landsat 8 TIRS images with 100 m resolution (resampled to 30 m) acquired on 10 April, 22 July, 3 October, and 22 December 2014 from the USGS (United States Geological Survey) website, which represented spring, summer, autumn, and winter, respectively, were utilized to derive LST to characterize the urban thermal environment. There was almost no cloud cover in the study area on the dates of acquisition. A mono-window algorithm (MWA) presented by Qin et al. [45] was used to map LST, with only three parameters required: emissivity, transmittance, and effective mean atmospheric temperature without the use of atmospheric profiles. The basic form can be written as: where T s is LST, T i is brightness temperature of TIRS bands i, and T a is the effective mean atmospheric temperature. The value of constants a and b are −67.366351 and 0.458606, respectively, and C and D are defined as follows: where ε and τ are land surface emissivity (LSE) of band i and atmospheric transmittance of band i, respectively. In this study, Landsat 8 TIRS band 10 takes the place of TM band 6 in MWA. The calculations for T i , ε, τ, and T a are detailed in previous studies [46][47][48][49][50]. Air temperature collected from 19 local meteorological stations located in the study area were used to validate the precision of the retrieved LST from Landsat images. The acquisition time of the two remote sensing images on 31 July 2014 was 10:21 a.m. However, measurements of the ground-truth LST were limited by the difficulty of finding a homogeneous region as large as the satellite pixel size. As a result, the average retrieved LST of the 9 × 9 pixels near the meteorological station was compared to the mean air temperature acquired on 10:00 a.m. and 11:00 a.m.

Characterizing the Urban Green Spaces
Landscape ecology, if not ecology in general, is largely founded on the notion that environmental patterns strongly influence ecological processes [51]. Landscape metrics are versatile and widely used to measure landscape patterns [40,43]. Six landscape metrics were chosen to quantitatively characterize the urban green spaces based on the following principles: (1) easily calculated; (2) scientific and practical importance [52,53]. The patch level metrics were patch area (PA), patch perimeter (PERIM), perimeter-area ratio (PARA), shape index (SHAPE), and class level metrics, including total area (TA). The landscape level metrics were: total area (TA) and number of patches (NP). All of these metrics and their descriptions are listed in Table 2.
As composition metrics, PA and TA are perhaps the most important and useful pieces of information regarding the landscape because they are used as the basis for many of the patch, class, and landscape indices [43]. PARA and SHAPE are shape metrics that quantify landscape configuration in terms of the complexity of patch shape at different levels. For example, a square patch is SHAPE = 1, and a higher SHAPE value indicates higher shape complexity. NP is an aggregation metric that refers to the tendency of patch types to be spatially aggregated. This property is also often referred to as landscape texture. These metrics characterize the composition and configuration information of urban green spaces. In this study, all the metrics were calculated by FRAGSTATS, which is a spatial pattern analysis program for quantifying the structure of landscapes.
For landscape level metrics, in order to build quantitative relationships between urban green spaces and LST at multiple resolutions, two spatial scales were selected: 1 km × 1 km and 0.5 km × 0.5 km. Fishnet in ArcGIS 10.2 (ESRI, Redlands, CA, USA) was used to create these grids, and there were 577 and 2206 sections for the scales of 1 km × 1 km and 0.5 km × 0.5 km, respectively. Each section is treated as a single landscape, and landscape metrics and LST were calculated for every grid. As a result, the spatial pattern of the landscape level metrics could be shown as theme maps to provide more distinct information compared to conventional tables. Table 2. Landscape metrics used in this study to describe urban green spaces.

Metrics and Abbreviation
Descriptions and Calculations patch area, PA equals the area (m 2 or hectare) of the patch, PA > 0, without limit patch perimeter, PERIM equals the perimeter (m) of the patch, PERIM > 0, without limit perimeter-area ratio, PARA equals the ratio of the patch perimeter (m) to area (m 2 ). PARA = PERIM/PA, PARA > 0, without limit shape index, SHAPE equals the ratio of the patch perimeter (m) to area (m 2 ). SHAPE = 0.25 PERIM/ √ PA, SHAPE ≥ 1 total area, TA equals the sum of the areas (m 2 ) of all patches of the corresponding patch type or the landscape, TA > 0, without limit number of patches, NP equals the number of patches of the corresponding patch type (class), NP ≥ 1, without limit In addition, at the patch and class level, NDVI (normalized difference vegetation index) was used to describe the amount of urban green spaces and growth status. High values of NDVI tend to correlate with more vegetation and good "condition" of vegetation. NDVI is defined as follows: where NIR is band 5 and R is band 4 for Landsat 8.

Statistical and Comparative Analysis
Existing research studies on relationships between urban green spaces and LST have used both qualitative and quantitative models [9,50]. In this study, all statistical analyses were conducted with SPSS 19.0. First, one-way analysis of variance (ANOVA) was conducted to examine whether the landscape metrics calculated from the three levels, the LST and NDVI of each patch, class, and landscape were significantly different and redundant. Then, bivariate correlations were performed to analyze the relationships between these metrics and corresponding LST. At the patch level, as the distribution of the area of the urban green space patches is not normal, Spearman's rho correlation coefficients were chosen to explore the relationships between the mean LST for each patch, landscape metrics, and NDVI on 31 July 2014. The area of a single TIRS image pixel was 900 m 2 (0.09 ha). In cases where the PA is not as large as the area depicted in the pixel, it is difficult to obtain accurate LST for each patch. For this reason, in our case, only the green space patches with areas larger than 0.36 ha (≥4 TIRS pixels) were selected to perform bivariate correlations, and these were classified into six classes depending on the patch size. At the landscape level, eight regression models were performed using the metrics and mean LST of each grid across the four seasons at different scales to identify the metrics most relevant to the LST.

Spatial Pattern of Urban Green Spaces
The spatial distribution of urban green spaces in Changchun is shown in Figure 2. The total area of green space is distributed across 115.84 km 2 with 8287 patches, and all are located within the 5th Ring Road. As shown in Table 3, Attached Green Spaces (AGS) patches were observed to have the smallest mean area, while Park Green Spaces (PGS) patches had the largest. However, AGS was observed to have the largest total area followed by Ecological Green Spaces (EGS), Road Green Spaces (RGS), PGS, Other Green Spaces (OGS), and Productive Plantation Green Spaces (PPGS). Table 3 provides detailed statistical information about the urban green spaces.

Spatial Pattern of LST
The LST distribution in Figure 3 shows the contrasts between surface temperature for different dates. The results showed that the average satellite-based LST was approximately 4.6 °C higher than the average air temperature, which was reasonable for the hot summer. Hence, the retrieved LST by the same algorithm in other seasons can also reflect the urban thermal environment. There were remarkable seasonal variations in LST. The average LST of the study area was 20.97 °C, 35.57 °C, 16.35 °C, and −14.56 °C in spring, summer, autumn, and winter, respectively. In spring, summer, and autumn, lower surface temperatures were found in urban green spaces. In contrast, in the winter, this phenomenon was not as apparent. The largest range of LST was found in the summer, representing the largest temperature contrast of all of the seasons. Detailed information of the LST in the study area for the four seasons is shown in Table 4.

Spatial Pattern of LST
The LST distribution in Figure 3 shows the contrasts between surface temperature for different dates. The results showed that the average satellite-based LST was approximately 4.6 • C higher than the average air temperature, which was reasonable for the hot summer. Hence, the retrieved LST by the same algorithm in other seasons can also reflect the urban thermal environment. There were remarkable seasonal variations in LST. The average LST of the study area was 20.97 • C, 35.57 • C, 16.35 • C, and −14.56 • C in spring, summer, autumn, and winter, respectively. In spring, summer, and autumn, lower surface temperatures were found in urban green spaces. In contrast, in the winter, this phenomenon was not as apparent. The largest range of LST was found in the summer, representing the largest temperature contrast of all of the seasons. Detailed information of the LST in the study area for the four seasons is shown in Table 4. Table 4. The statistical information of the land surface temperature (LST) among the four seasons.

Date
Max

Patch Level
As an indicator of the status of vegetation, Table 5 shows that the mean NDVI of each patch was significantly negatively correlated to LST among all the classes, indicating that higher NDVI values can lead to lower LST. When the area of green space patch is between 1.44 and 5.76 ha, the negative relationship between LST and PA is significant and much stronger than for the other classes, indicating that in this area range increasing the PA can lower LST. While in other area ranges, it appears that the PA has little influence on its own LST. For other landscape metrics, the relationships between LST and PERIM, PARA, and SHAPE were very weak. There are two exceptions. When the PA is larger than 15.2 ha, the negative relationships between LST and PERIM, and SHAPE were

Patch Level
As an indicator of the status of vegetation, Table 5 shows that the mean NDVI of each patch was significantly negatively correlated to LST among all the classes, indicating that higher NDVI values can lead to lower LST. When the area of green space patch is between 1.44 and 5.76 ha, the negative relationship between LST and PA is significant and much stronger than for the other classes, indicating that in this area range increasing the PA can lower LST. While in other area ranges, it appears that the PA has little influence on its own LST. For other landscape metrics, the relationships between LST and PERIM, PARA, and SHAPE were very weak. There are two exceptions. When the PA is larger than 15.2 ha, the negative relationships between LST and PERIM, and SHAPE were significant and stronger. It appears that the SHAPE metric begins to have influence on the mean LST of green space patches in this area range. These relationships suggest that compared to NDVI, landscape metrics including PERIM, PARA, and SHAPE have little influence on the mean LST of green space patches. When the PA of one patch is less than 1.44 ha or larger than 5.76 ha, it has little effect on its own mean LST. Table 5. Spearman's rho correlation coefficients between LST and patch-level metrics, normalized differential vegetation index (NDVI) on 31 July 2014. A summary of the stepwise regression model is shown in Table 6. The results showed that the standardized coefficients of variables differed greatly and all variables except for PA were included in the best fit model, indicating that the PA of each patch had little influence on its own mean LST. Compared to other patch-level metrics, NDVI is the most important variable in determining the LST. Different from other variables, a positive relationship was detected between PARA and LST.

Class Level
The average LST and standard deviation of each green type were lower than the mean LST of the city as a whole in all dates except for 22 December, as shown in Tables 4 and 7. The NDVI for each green type was greater than the average values of the overall area with 0.12, 0.42, 0.25, and 0.05 for 10 April, 31 July, 3 October, and 22 December, respectively. For seasonal variations, the correlation coefficients between the mean LST and NDVI of the six green spaces classes were −0.42, −0.94 (0.01 level, one-tailed), −0.93 (0.01 level, one-tailed), and −0.74 (0.05 level, one-tailed) for dates 10 April, 31 July, 3 October, and 22 December, respectively. PGS had the lowest mean LST and highest NDVI across all seasons except winter, while AGS and RGS had a high mean LST and a low NDVI relative to that of PGS.

Landscape Level
The temperature difference (TD) and NDVI difference between the whole study area and urban green spaces across the four dates are shown in Figure 4. In the winter, there was almost no TD between urban green spaces and the entire study area. The NDVI values of green spaces were lower than those of the study area across all seasons with the greatest difference in summer and the least in winter.

Landscape Level
The temperature difference (TD) and NDVI difference between the whole study area and urban green spaces across the four dates are shown in Figure 4. In the winter, there was almost no TD between urban green spaces and the entire study area. The NDVI values of green spaces were lower than those of the study area across all seasons with the greatest difference in summer and the least in winter. The spatial distributions of landscape-level metrics TA and NP are shown in Figure 5. Higher vegetation cover is located in the northeast and southeast part of study area where there are several big parks or lots of AGS patches in high density. The number of patches (NP) is high within the 4th Ring Road because there are many small AGS patches located in the oldest parts of the city. The spatial distributions of the two metrics suggests there is no obvious relationship between them. Further quantitative regression models were established to explore the relationships between the landscape-level metrics and LST at 0.5 km × 0.5 km and 1 km × 1 km scales. Two-dimensional scatter plots are shown in Figures 6 and 7. Negative linear relationships between LST and TA were found on all dates except for 22 December at the two different spatial scales. For season variations, negative relationships were strongest in summer followed by autumn and spring, and the relationships in winter were not significant. For spatial scale differences, compared to 0.5 km × 0.5 km, the negative relationships were stronger for the 1 km × 1 km scales across all the seasons. On 31 July at the 0.5 km × 0.5 km scale, a 1 ha increase in TA can result in approximately a 0.34 °C decrease in temperature. In other words, a 10% increase in TA can lead to a 0.85 °C decrease in temperature The spatial distributions of landscape-level metrics TA and NP are shown in Figure 5. Higher vegetation cover is located in the northeast and southeast part of study area where there are several big parks or lots of AGS patches in high density. The number of patches (NP) is high within the 4th Ring Road because there are many small AGS patches located in the oldest parts of the city. The spatial distributions of the two metrics suggests there is no obvious relationship between them.

Landscape Level
The temperature difference (TD) and NDVI difference between the whole study area and urban green spaces across the four dates are shown in Figure 4. In the winter, there was almost no TD between urban green spaces and the entire study area. The NDVI values of green spaces were lower than those of the study area across all seasons with the greatest difference in summer and the least in winter. The spatial distributions of landscape-level metrics TA and NP are shown in Figure 5. Higher vegetation cover is located in the northeast and southeast part of study area where there are several big parks or lots of AGS patches in high density. The number of patches (NP) is high within the 4th Ring Road because there are many small AGS patches located in the oldest parts of the city. The spatial distributions of the two metrics suggests there is no obvious relationship between them. Further quantitative regression models were established to explore the relationships between the landscape-level metrics and LST at 0.5 km × 0.5 km and 1 km × 1 km scales. Two-dimensional scatter plots are shown in Figures 6 and 7. Negative linear relationships between LST and TA were found on all dates except for 22 December at the two different spatial scales. For season variations, negative relationships were strongest in summer followed by autumn and spring, and the relationships in winter were not significant. For spatial scale differences, compared to 0.5 km × 0.5 km, the negative relationships were stronger for the 1 km × 1 km scales across all the seasons. On 31 July at the 0.5 km × 0.5 km scale, a 1 ha increase in TA can result in approximately a 0.34 °C decrease in temperature. In other words, a 10% increase in TA can lead to a 0.85 °C decrease in temperature Further quantitative regression models were established to explore the relationships between the landscape-level metrics and LST at 0.5 km × 0.5 km and 1 km × 1 km scales. Two-dimensional scatter plots are shown in Figures 6 and 7. Negative linear relationships between LST and TA were found on all dates except for 22 December at the two different spatial scales. For season variations, negative relationships were strongest in summer followed by autumn and spring, and the relationships in winter were not significant. For spatial scale differences, compared to 0.5 km × 0.5 km, the negative relationships were stronger for the 1 km × 1 km scales across all the seasons. On 31 July at the 0.5 km × 0.5 km scale, a 1 ha increase in TA can result in approximately a 0.34 • C decrease in temperature. In other words, a 10% increase in TA can lead to a 0.85 • C decrease in temperature on this spatial scale. While on the same date at a 1 km × 1 km scale, a 10% increase in TA can lead to a 1.2 • C decrease in temperature, according to the statistical analysis. In other words, with an increase in TA, urban green spaces can result in a temperature decrease, but the cooling effect varies greatly at different seasons and spatial scales. on this spatial scale. While on the same date at a 1 km × 1 km scale, a 10% increase in TA can lead to a 1.2 °C decrease in temperature, according to the statistical analysis. In other words, with an increase in TA, urban green spaces can result in a temperature decrease, but the cooling effect varies greatly at different seasons and spatial scales.   on this spatial scale. While on the same date at a 1 km × 1 km scale, a 10% increase in TA can lead to a 1.2 °C decrease in temperature, according to the statistical analysis. In other words, with an increase in TA, urban green spaces can result in a temperature decrease, but the cooling effect varies greatly at different seasons and spatial scales.   Three-dimensional images of TA (ha), NP, and LST across different seasons at different spatial scales were generated to explore the effect of configuration of urban green spaces on their cooling intensity. For the same area of urban green spaces, it appears that their spatial configurations can change significantly. For example, the number of patches and their locations play an important role in decreasing LST. For urban green space planning, the configuration cannot be ignored. In this study, the NP in particular spatial scales was chosen as the indicator to characterize the information for the configuration. Figure 8 shows the 3-D images at the 0.5 km × 0.5 km scale. On 10 April, when the TA is the same, with the increase of NP, the LST has a negative trend. On 31 July, when the TA is less than 5 ha, the NP has little influence on LST. When TA is between 5 and 15 ha, an increase of NP can lead to lower LST. When TA is larger than 15 ha at this scale, NP is no longer the key factor in determining LST. On 3 October, when TA is less than 10 ha, NP has a negative relationship with LST. However, when TA is larger than 10 ha, the effect of NP on LST is not significant. On 22 December, because the urban green spaces have little influence on LST, their configurations also have little impact on LST. Figure 9 shows the 3-D images at the 1 km × 1 km scale. On 10 April, with an increase of NP, the LST tends to have lower values. On 31 July, when the proportion of TA is less than 20%, regardless of NP, the LST is relatively high. When the proportion of TA is between 20% and 40%, at the beginning, with an increase of NP, the LST increases simultaneously. However, when NP is larger than 40, the LST decreases. When TA is larger than 40%, less NP can have a better cooling effect. On 3 October, when the percentage of TA is less than 20%, if NP is smaller than 50, the NP has little influence on LST. However, NP can result in a lower LST when it is larger than 50. A very interesting phenomenon is that with the increase of TA, the threshold value of NP can lead to lower LST decreases. For example, when the TA is between 30 and 40 ha, the threshold of NP is approximately 30. In contrast, when the TA is between 50 and 60 ha, the threshold of NP is approximately 10. However, when the NP is larger than the threshold value, the increase of NP has little effect on LST. These results indicate that when the TA is larger in this spatial scale, less NP may generate a better cooling effect. Compared to the 0.5 km × 0.5 km scale, on 22 December, with the increase of NP, the LST had an increasing trend at 1 km × 1 km scale. Three-dimensional images of TA (ha), NP, and LST across different seasons at different spatial scales were generated to explore the effect of configuration of urban green spaces on their cooling intensity. For the same area of urban green spaces, it appears that their spatial configurations can change significantly. For example, the number of patches and their locations play an important role in decreasing LST. For urban green space planning, the configuration cannot be ignored. In this study, the NP in particular spatial scales was chosen as the indicator to characterize the information for the configuration. Figure 8 shows the 3-D images at the 0.5 km × 0.5 km scale. On 10 April, when the TA is the same, with the increase of NP, the LST has a negative trend. On 31 July, when the TA is less than 5 ha, the NP has little influence on LST. When TA is between 5 and 15 ha, an increase of NP can lead to lower LST. When TA is larger than 15 ha at this scale, NP is no longer the key factor in determining LST. On 3 October, when TA is less than 10 ha, NP has a negative relationship with LST. However, when TA is larger than 10 ha, the effect of NP on LST is not significant. On 22 December, because the urban green spaces have little influence on LST, their configurations also have little impact on LST. Figure 9 shows the 3-D images at the 1 km × 1 km scale. On 10 April, with an increase of NP, the LST tends to have lower values. On 31 July, when the proportion of TA is less than 20%, regardless of NP, the LST is relatively high. When the proportion of TA is between 20% and 40%, at the beginning, with an increase of NP, the LST increases simultaneously. However, when NP is larger than 40, the LST decreases. When TA is larger than 40%, less NP can have a better cooling effect. On 3 October, when the percentage of TA is less than 20%, if NP is smaller than 50, the NP has little influence on LST. However, NP can result in a lower LST when it is larger than 50. A very interesting phenomenon is that with the increase of TA, the threshold value of NP can lead to lower LST decreases. For example, when the TA is between 30 and 40 ha, the threshold of NP is approximately 30. In contrast, when the TA is between 50 and 60 ha, the threshold of NP is approximately 10. However, when the NP is larger than the threshold value, the increase of NP has little effect on LST. These results indicate that when the TA is larger in this spatial scale, less NP may generate a better cooling effect. Compared to the 0.5 km × 0.5 km scale, on 22 December, with the increase of NP, the LST had an increasing trend at 1 km × 1 km scale.

Discussion
Representative landscape metrics at patch, class, and landscape level and quantitative statistical analysis were utilized to explore the effect of urban green spaces on the urban thermal environment. The results demonstrated that not only the composition of urban green spaces but also their configuration can significantly affect the distribution of the LST. The effect, however, varies greatly across different seasons depending considerably on the size and shape of urban green spaces, the status of vegetation, the analysis spatial scales, and seasonal radiation conditions. The results were consistent with many other previous studies, which demonstrated negative correlations between LST and urban green spaces [3,19,54].

Implications for Methodology
First, comparative approaches with an additional cold zone case study were utilized in this study to investigate the relationship between urban green spaces and the urban thermal environment. Comparative approaches included a selection of landscape metrics at three levels, patch, class, and landscape level. The four dates representing spring, summer, autumn, and winter were chosen to explore seasonal variations; two spatial scales, 0.5 km × 0.5 km and 1 km × 1 km, were selected to study the relationship changes across scales. These helped generalize a more thorough scientific understanding of the effects of urban green spaces on LST.
Second, we propose that spatialization is an important part of landscape metrics. Traditionally, the result of landscape metrics calculated by metrics software, such as FRAGSTATS, is often shown in tables. In this study, with the help of ArcGIS 10.2, 2206 and 577 grids were created in the study area at the scale of 0.5 km × 0.5 km and 1 km × 1 km, respectively. Each grid was considered as a whole landscape, and landscape-level metrics were calculated for each grid. As a result, the spatial information of the metrics were presented as maps in order to show the distribution of metrics in a more explicit way.

Theoretical Implications
In this study, landscape metrics and NDVI were used to characterize the composition features of urban green spaces.

Discussion
Representative landscape metrics at patch, class, and landscape level and quantitative statistical analysis were utilized to explore the effect of urban green spaces on the urban thermal environment. The results demonstrated that not only the composition of urban green spaces but also their configuration can significantly affect the distribution of the LST. The effect, however, varies greatly across different seasons depending considerably on the size and shape of urban green spaces, the status of vegetation, the analysis spatial scales, and seasonal radiation conditions. The results were consistent with many other previous studies, which demonstrated negative correlations between LST and urban green spaces [3,19,54].

Implications for Methodology
First, comparative approaches with an additional cold zone case study were utilized in this study to investigate the relationship between urban green spaces and the urban thermal environment. Comparative approaches included a selection of landscape metrics at three levels, patch, class, and landscape level. The four dates representing spring, summer, autumn, and winter were chosen to explore seasonal variations; two spatial scales, 0.5 km × 0.5 km and 1 km × 1 km, were selected to study the relationship changes across scales. These helped generalize a more thorough scientific understanding of the effects of urban green spaces on LST.
Second, we propose that spatialization is an important part of landscape metrics. Traditionally, the result of landscape metrics calculated by metrics software, such as FRAGSTATS, is often shown in tables. In this study, with the help of ArcGIS 10.2, 2206 and 577 grids were created in the study area at the scale of 0.5 km × 0.5 km and 1 km × 1 km, respectively. Each grid was considered as a whole landscape, and landscape-level metrics were calculated for each grid. As a result, the spatial information of the metrics were presented as maps in order to show the distribution of metrics in a more explicit way.

Theoretical Implications
In this study, landscape metrics and NDVI were used to characterize the composition features of urban green spaces.
The negative relationship between urban green space areas (patch area, PA) and LST observed in this study is in agreement with previous studies [55][56][57]. A higher percentage of cover of green spaces leads to lower LST, and its negative effects on LST is more significant in summer than for the other three seasons. This may be due to summer having the most heterogeneous thermal environment and higher evapotranspiration of vegetation. Compared to the area of urban green patches, the negative effect of PERIM on LST is much weaker, indicating that it may not be an ideal metric to present information of green spaces for this study. We found that shape metrics PARA and SHAPE, compared to area metrics (PA, TA), were not very valuable in all seasons for assessing their relationship with LST. Moreover, the results suggested that a significant cooling effect occurrs when the TA is beyond a certain threshold. For example, in summer at the 0.5 km × 0.5 km scale, the value is approximately 5 ha. Therefore, based on the results, increasing the area of urban green spaces can effectively mitigate the UHI, and this was consistent with previous studies [38,58].
NDVI is typically used as the indicator of vegetation abundance to estimate the LST-vegetation relationship [21,59]. The results indicated that high values of NDVI for urban green spaces led to significantly lower temperatures for all seasons except winter at class level. The negative relationship between LST and NDVI was much stronger than that of LST and the area of urban green spaces. High values of NDVI were usually related to a good grow status of vegetation or more green space density. This is a clear indication that the status of vegetation (e.g., tree species, vegetation cover percentage) plays a vital role in cooling the urban thermal environment. In winter, the urban green spaces have little effect on the urban thermal environment, confirming that the status of vegetation does matter for cooling intensity. Therefore, in an urban environment without enough public space for green vegetation, increasing the NDVI (for example, vegetation canopy cover) is still an effective way to decrease LST.
For different types of urban green spaces, PGS has a much lower temperature than the other types, while AGS and RGS usually have relative high temperatures. This may be due to several reasons: (1) the sizes of PGS sites were larger than other green space types, and there were more trees in parks; (2) the AGS patches were usually small and showed discrete distributions; and (3) the RGS sites were isolated and there were anthropogenic heat sources beside them, such as vehicles and buildings. As a result, by only considering the cooling intensity, if there is enough space, parks are a better choice than AGS and RGS for urban green infrastructure.
The whole study area was divided into small grids with a certain size, therefore the landscape-level metric NP rather than PD (patch density) was used to assess the spatial pattern of urban green spaces on LST. Our results showed that although the effect of the configuration of green spaces was not as strong as compositional features, it did play an important role in urban green space planning for regulating the micro-climate. For example, at the 1 km × 1 km scale, if a certain area of green space was to be arranged, how many patches should be considered? Our results showed that if the total area of green spaces was not large enough (approximately 20 ha in our case), the configuration has little influence. In summer, when the total area is larger than the threshold, for example, when the TA is between 20 and 30 ha in this study, with an increase in the number of patches, the LST tends to have few variations and even increases. However, when the patch number is larger than the threshold, the LST of the 1 km × 1 km grid will decrease significantly. Although when the TA is large enough on this scale, the NP has little effect on LST. As a result, the effects of NP on LST depends on the total area of urban green spaces and the season. Therefore, it still remains a challenge to assess how the NP can decrease LST in the most effective way. In addition, other landscape-level metrics, such as splitting index or connectance index could also be employed to investigate the effect of configuration in future studies. The splitting index is usually used to characterize landscape fragmentation in a geometric perspective and increases as the landscape is increasingly subdivided into smaller patches [60]. The connectance index is used as the proportion of functional joinings among all patches, where each pair of patches is either connected or not based on some criterion [43].

Limitations
This study has several limitations. First, only Landsat 8 TIR data were used to retrieve LST to characterize the spatial patterns of the urban thermal environment, thus the satellite overpass time was not synchronous with that of sparsely distributed meteorological stations. Due to the 30-m spatial resolution of TIRS data, the LST of some small green patches was influenced by surrounding urban areas. In addition, only one date was chosen to represent a season because there were insufficient daily or monthly images to meet the demand for retrieving LST. In addition, it is well known that different atmospheric conditions can affect the accuracy of the retrieved LST [26]. However, these limitations were related to the difficulties in retrieving LST.
The second limitation is that bodies of water also play a very important role in cooling the LST [37,61]. For example, if bodies of water extensively exist in one analysis grid, the cooling effect of urban green spaces may be overestimated. However, the cooling effect of water was not considered in our study because we only sought to explore the effect of urban green spaces on thermal environments at the landscape level.
Finally, more attention is needed on the use of synchronous data from in-situ observations such as air temperature/humidity combined with high temporal resolution and long time span satellite data to have a better understanding of the cooling effect of urban green spaces produced by future studies. Researchers should interpret the results in comparison to previous studies and the working hypotheses. The findings and their implications should be further discussed in the broadest context possible, particularly to highlight future research directions.

Conclusions
Using Changchun, a cold temperate zone city, as a case study, this study investigated the effects of urban green spaces on urban thermal environments and their seasonal variations. With the support of multi-temporal Landsat OLI/TIRS images, high spatial resolution images GF-1, field data, and landscape ecology and comparative approaches were used to quantitatively examine the relationship between urban green spaces and LST. Our results showed the following: (1) Urban green spaces have a significant cooling effect among all seasons except for winter, but the effects vary considerably among different seasons and green space types. The difference of LST between the urban green spaces and urban areas was large in the summer and small in the winter. The maximum value of the differences was 1.27 • C on 31 July. The lowest LST was found in PGS, while AGS and RGS have the highest LST among different green patches. (2) Compared to shape metrics such as PARA and SHAPE, the negative relationships between LST and the area and NDVI of urban green spaces were more significant. We found that every 10% increase in urban green spaces explained a 0.34 • C decrease in LST in summer at the 0.5 km × 0.5 km scale. Rather than increasing the area of green patches, increasing the NDVI (more woodland area and denser vegetation) is a better and more practicable approach for urban green planning. (3) Both the compositional features of urban green spaces and configuration can significantly affect LST. At certain spatial scales, for example, 1 km × 1 km, given a fixed area of urban green spaces, the number of green patches can increase or decrease the values of LST depending on if the number is larger than a threshold or not, and this threshold changes for each given area. For example, in our case, when the proportion of TA is between 20% and 40% in 100 ha, initially, the LST increases simultaneously with an increase in NP. When NP is larger than 40, however, the LST decreases.
This study has the potential to contribute to our understanding of how urban green spaces can mitigate urban heat island effects so that environmental problems may be resolved. The outcome of this study suggested that the process of cooling effects of urban green spaces differed among different seasons and green types, and it may provide insights for city planners and decision makers to forge sustainable urban strategies by considering the composition and configuration features of urban green spaces at the same time.