Zhengzhou (34.16–34.58 N, 112.42–114.14 E) is located in the north temperate zone and has a temperate continental monsoon climate with four distinct seasons. The annual average temperature is 14.2°C, and the annual rainfall is 640.9 mm, which is concentrated in summer. As the capital city of Henan province, Zhengzhou is also a transportation hub in China. By 2021, the resident population of Zhengzhou was about 12.7 million, and the urbanization rate was 78.4%. The built-up area of the central urban area is 709.69 km². Here, the central city of Zhengzhou is taken as the study area, with a total area of 1019.5 km². A total of 36 parks within the urban center were selected (Figure 2, attached Table 1).

Landsat 8 remote sensing images taken at 11:01:13 on July 7, 2019, were utilized (http://www.gscloud.cn/), with good imaging quality and clear ground features and no clouds and band stripe. The remote sensing imager is equipped with Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS), including 11 bands, among which two thermal infrared bands can be used to retrieve land surface temperature (LST) data. The first few days of the study period were sunny and windless, and the weather conditions were good, excluding the influence of factors such as rainfall on the results.

In order to reduce or eliminate the influence of atmospheric environmental conditions on remote sensing image quality, radiometric calibration and atmospheric correction were carried out on multi-spectral bands. In addition, the corrected Landsat 8 remote sensing data were classified by the support vector machine method. According to the research needs, three categories of green space, water body, and impervious water surface were obtained. The classification accuracy was 99%, and the verified Kappa coefficient was 0.98, which met the application requirements. At the same time, parks in the city are extracted using vector boundary data for clipping (http://www.amap.com). Data processing was carried out in ENVI 5.3 and ArcGIS 10.8.

The PCS from the park to its vicinity was gradually mitigated by the urban surface heat environment (Jansson et al., 2006). During daytime, the impervious surface has a fast heating rate, resulting in the air temperature and LST being similar in the park, where the green spaces have relatively small thermal inertia and a fast heating rate, which causes the cooling effect in the park to be less pronounced (Yan et al., 2018). Here, we used the LST retrieved from remote sensing data instead of the air temperature. Generally, air temperature and LST are mutually influential, and several studies have found evidence of a close statistical relationship and similar spatial patterns between them (Schwarz et al., 2012). LST has been applied as a proxy to study air temperature. Our study requires a large sample of parks to explore solutions to enhance PCS, and in this case, it is more convenient and reasonable to use the LST obtained from the satellite infrared band.

At present, LST inversion algorithms mainly include the atmospheric correction method, single-window algorithm, and split window algorithm (Yu et al., 2014; Qiu and Jia, 2020; Liao et al., 2021). Referring to other research studies and according to the characteristics of the study area (Cao et al., 2010; Li et al., 2012; Sun et al., 2012; Li et al., 2020), we used the atmospheric correction method to invert land surface temperature. The atmospheric correction method, also known as the radiative transfer equation (RTE), is a traditional algorithm based on the atmospheric radiation transfer model. The principle is to remove the errors caused by atmospheric heat radiation to the surface from the total heat radiation observed by satellite sensors (Li et al., 2012; Barsi et al., 2014). The thermal radiation intensity is then converted to the corresponding surface temperature. The formula is as follows (Wang et al., 2018; Xue et al., 2019): T S = K 2 ln K 1 B T S + 1 . (1) B T S = L λ − L ↑ − τ 1 − ε L ↓ τ ε . (2) Here, T S is the surface temperature, B T S is the brightness value of blackbody radiation, and K 1 and K 2 are the coefficients. In TIRS 10 band of Landsat 8, K 1 =774.885 W ∙ m − 2 ∙ s r − 1 ∙ μ m − 1 and K 2 =1321.079 W ∙ m − 2 ∙ s r − 1 ∙ μ m − 1 . L ↑ and L ↓ are the upward and downward radiation intensity to the atmosphere, respectively, and τ is the atmospheric transmittance. The three parameters can be obtained from NASA’s website (http://atmcorr.gsfc.nasa.gov). ε is the surface emissivity, and L λ is the thermal infrared brightness value received by the sensor. ε can be calculated by the normalized vegetation index (NDVI) using the following methods (Sobrino et al., 2004): ε = 0.004 P V + 0.986 . (3) P V = N D V I − N D V I s o i l N D V I v e g − N D V I s o i l   0.05 ≪ N D V I ≪ 0.75 . (4) Here, P V is the vegetation coverage and N D V I v e g and N D V I s o i l are NDVI of vegetation and bare land, respectively. In addition, when NDVI exceeds 0.75, P V equals 1; when NDVI is less than 0.05, P V equals 0.

In order to quantify the PCS, we selected the spillover distance of the park cooling effect (PCS_D), spillover strength of the park cooling effect (PCS_S), and spillover rate of the park cooling effect (PCS_R) to describe it. Previous studies have shown that PCS_D is concentrated within 1000 m (Chang et al., 2007; Xie and Li, 2020; Zhu et al., 2021). To guarantee that PCS_D is within the buffer distance, we set up multi-stage buffer zones in 30-m units along the park boundary within 1500 m outside the park and statistically obtained the average LST within each buffer zone T₁, T₂, and T₃. The piecewise function was used for fitting analysis of these values, and the LST curve model around the park was established (Figure 3).

L, PCS_D, is the horizontal distance between the park boundary and the first inflection point of the LST curve (Qiu and Jia, 2020), which is located at the position where the slope of the LST curve changes sharply or tends to be stable (Figure 3) (Sun et al., 2012; Du et al., 2016). PCS_S is the difference between the LST at the cooling boundary of the park and the average LST inside the park, and the calculation formula is as follows: P C S S = T L − T m e a n . (5) Here, T L is the LST at which the park reaches its maximum PCS_D, and T m e a n is the average LST inside the park. PCS_R is the LST change rate within unit distance, and the calculation formula is as follows: P C S R = P C S s P C S D . (6)

The fitting of the piecewise function and the determination of the inflection point were completed in R 4.1.3’s segmented package.

In order to analyze the relationship between the landscape pattern inside the park and PCS, this study selected 11 indicators (Table 1) for analysis from the overall characteristics, internal composition, and landscape shape of the park. The index calculation was completed in FRAGSTATS 4.2 and ArcGIS 10.8, and the significance and calculation method of each index can be referred to Wu (2013).

PCS is a non-linear process. In the context of scarce urban land resources, it is essential to quantify the threshold of park cooling efficiency and determine the optimal park size. The maximization method was adopted to standardize the park area and PCS index, and the logarithmic function was used for fitting. The formula is as follows: y = a ∙ l n x + b . (7) y ′ = a ∙ x − 1 . (8) Here, y is PCS (PCS_D, PCS_S, and PCS_R), x is the landscape index related to park area (AREA; AREA_g), and y ′ is the cooling efficiency, which is optimal when the slope is equal to 1 (Peng et al., 2020). Before the threshold, a small increase in the park area will lead to a considerable increase in PCS, while when the park area exceeds the threshold, PCS cannot be significantly enhanced by increasing the park area (Zhu et al., 2021; Yao et al., 2022). So, the corresponding park area when the slope of the fitting curve is 1 can be identified as the threshold.

Based on formula (1) ∼ (4), the average LST in the study area is 37.8°C, which is consistent with field observation for noon temperature in summer. The LST in the city center is apparently higher than that in the suburbs (Figure 4A). The average LST in the First Ring Road is 1.4°C higher than that outside the Fourth Ring Road, indicating an obvious UHI effect (Figure 4B). Due to the different characteristics of the underlying surface, the fluctuation of LST in UHI is significant. The main reason is that some urban parks, rivers, and lakes have low LST, forming obvious low-temperature patches, such as Longhu Park and Xiliuhu Park, which play the role of a “cold island” in the city.

The spatial superposition analysis of park patches and LST showed that the average LST in 36 parks is 35.5°C, which is 2.3°C lower than the average LST in Zhengzhou, indicating that there is a clear cooling effect in urban parks. However, the average LST in different parks varies significantly (Figure 5). There are 28 parks having LST lower than the average LST in Zhengzhou, among which Park 2 and Park 3 have the lowest LST (33.9°C). The LST of Park 34 is the highest, reaching 40.7°C. In addition, the LST values of Park 31 and Park 36 are both higher than 38.8°C, exceeding the average LST of Zhengzhou.

The average LST of the multi-stage buffer zone in each park was extracted, and the LST varied with distance in four types (Figure 6). The changes of LST in the buffer zone of 33 parks are similar to those given in Figures 6A–C, such as Park 3, Park 18, and Park 26. The LST–distance fitting curve of such parks shows a trend of increasing first and then a flat trend, which indicates that the LST decreases significantly within a certain distance around them. PCS gradually weakens as the distance increases. When the distance exceeds its PCS_D, LST change is no longer obvious or fluctuates slightly due to the change in underlying surface properties (Li et al., 2013; Peng et al., 2020). On the contrary, in three parks (Park 34, Park 35, and Park 36), the LST in the buffer zone is lower than that in the park (Figure 6D), indicating that the park has a “thermal effect” spillover. We used piecewise function to perform fitting analysis on temperature values of multi-stage buffer and obtain PCS_D. The results showed that the fitting effect is good and all R² values are greater than 0.9.

Figure 7 shows the impact of the park on the surrounding environment. The average PCS_S of 33 parks with PCS in the study area is 1.8°C, where the PCS_S values of Park 1 and Park 2 are the highest, with both exceeding 4°C. PCS_S values of Park 32 and Park 33 are the lowest, being only 0.04°C and 0.05°C, respectively. Meanwhile, the difference in PCS_D is significant, with an average of 179 m (Table A2). The number of parks with PCS_D between 100 and 200 m forms the majority, accounting for 39.4% of the total number of samples. There are only six parks with PCS_D more than 300 m, accounting for 18.2%. PCS_D values in Park 1 and Park 4 are larger, at 370 m and 345 m, respectively. The smallest PCS_D is of Park 32, which is only 31 m. Overall, parks have the greatest cooling effect on the surrounding area within 200 m. The average PCS_R of 33 parks with PCS is 0.96°C/100 m, which means that for every 100-m increase in distance from the park, the LST will increase by about 0.96°C. The PCS_R of Park 27 and Park 22 exceeded 1.70°C/100 m. The smaller PCS_R is of Park 33, with a cooling intensity of only 0.13°C per 100 m.

We analyzed the relationship between the 11 landscape characteristic indexes of the parks and PCS (PCS_D, PCS_S, and PCS_R) (Figure 8). The results showed that park landscape features have a strong influence on PCS. PARA, ED, ED_g, and PD were inversely correlated with PCS_D. The correlation coefficients were −0.81, −0.76, −0.72, and −0.57, respectively, which indicates that the increase of PARA and ED was not conducive to the generation of larger PCS_D. Qiu and Jia (2020) pointed out that PD is the most important negative indicator of the cooling effect. In this study, the increase in PD also has a negative effect on PCS_D, but the degree of influence is not as strong as that of the first three.

AREA, PERIM, AREA_g, LSI_i, LSI, LSI_g, and P_w were positively correlated with PCS_D. The correlation coefficients were 0.86, 0.85, 0.85, 0.79, 0.6, 0.59, and 0.55, respectively. Specifically, the role of AREA, PERIM, AREA_g, and LSI_i in promoting PCS_D is obvious, which is consistent with the research results of Cao et al. (2010) and Yao et al. (2022). Therefore, PCS_D can be significantly expanded by increasing AREA, AREA_g, and changing LSI_i, while the cooling benefits from changing LSI and LSI_g or increasing the proportion of water are less obvious.

In terms of PCS_S, ED and PARA had the most significant negative effect on it, with correlation coefficients of −0.74, followed by ED_g (−0.62) and PD (−0.56). Comparing the results of the PCS_D’s study, we found that the impacts of the four landscape indicators (PARA, ED, ED_g, and PD) on PCS_D and PCS_S are similar. The difference is that the correlation between these four indicators and PCS_S is not as good as that of PCS_D. This means that PARA, ED, ED_g, and PD have significant negative effects on PCS. How to alleviate the negative effects of these four aspects is crucial to improve the cooling effect of urban parks.

AREA_g, AREA, PERIM, LSI_i, and P_w were positively correlated with PCS_S, with correlation coefficients of 0.77, 0.76, 0.67, 0.6, and 0.57, respectively. The results demonstrated that AREA_g is the primary influencing factor and the degree of influence on PCS_S exceeds that of AREA and PERIM. LSI_i had a weaker effect on PCS_S than on PCS_D. P_w has a similar effect on both. In addition, LSI and LSI_g had no significant correlation with PCS_S.

PCS_R was not significantly correlated with landscape elements such as AREA and LSI but only had a positive correlation with PCS_S.

We used non-linear functions to fit the correlation between the area element of the park (AREA and AREA_g) and its PCS feature (PCS_D and PCS_S) (Figure 9). Formulas 7 and 8 were used to calculate the area threshold when the PCS efficiency of the park was at its highest. The results manifested that when the efficiency of the PCS_D is the highest, the best AREA was 8.08 ha (Figure 9A) and the best AREA_g was 5.71 ha (Figure 9B). The AREA with the best PCSS efficiency was 6.80 ha (Figure 9C), and the AREA_g was 5.08 ha (Figure 9D). In general, the park area to achieve high PCS efficiency in Zhengzhou city should be 6–8 ha, while the internal green space area should be controlled within 5–6 ha.

Urban parks are considered an effective way to mitigate the UHI (Chiesura, 2004; Ayala-Azcárraga et al., 2019), which is further confirmed by this study. The average LST of 36 parks in the study area is 2.3°C lower than that of Zhengzhou city, and the LST of the 33 parks is lower than that of the surrounding environment, forming an obvious “urban cold island.” However, the LST difference between different parks and the surrounding environment is not consistent, and the PCS_S is between 0.04°C and 4.61°C. The LST of only three parks was higher than that of the surrounding environment, forming a “heat island.”

The LST is closely related to the landscape characteristics within the park, such as the size of the green spaces and water bodies. Whether the spatial configuration of the patches facilitates ventilation and heat exchange entails further investigation. Existing studies have also demonstrated that the patch inside the larger ED park is severely fragmented, and it is arduous to form a relatively stable microclimate (Cheung et al., 2021; Wu et al., 2021). The larger the PARA, the stronger the heat exchange capacity is between the park and the surrounding environment (Yu et al., 2020), which is not conducive to the formation of low-temperature areas. Furthermore, the impervious water surface also seriously affects the change of the LST inside the park. As an example, in our study, the proportion of impervious water surface in the three parks with UHI effect is over 65%, and the mean LSI is 1.82, which is far lower than the mean 2.70 of 33 parks with PCS. This indicated that concentrated and continuous hardened surfaces increase park LST rather than decreasing it.

In order to explore the causes of thermal effects in parks (Park 34–36), we calculated the mean values of indicators for two types of parks and explored the relationship between park landscape indicators and park LST (Tables 2, 3). From the table, we can find that AREA_g, AREA, and LSI_i can effectively reduce LST, among which AREA_g is the most paramount influencing factor. However, these indicators of Parks 34–36 are much smaller than those of other parks. These three parks account for very little vegetation (34.8%), which means that they lack a cooling source and cannot efficiently reduce LST inside the park. Also, the mean values of PARA, PD, and ED in Parks 34–36 are higher than in other parks, and these three characteristics can cause the park temperature to increase. In this case, Parks 34–36 presented a situation where the park LST was higher than that of the surrounding environment.

PCS_D and PCS_S are major aspects of the cooling effect of the park. The larger PCS_D suggests that more individuals can enjoy the cooling benefits of parks, while the improvement of PCS_S means that the cooling effect directly felt by the human body is more obvious, both of which are very significant for improving the happiness of urban residents. We found the average PCS_S in the study area is 1.8°C, and the maximum can exceed 4°C. PCS_D can reach 370 m (Park 1) and has the most obvious cooling effect on the surrounding 200 m, and the average PCS_R is 0.957°C/100 m. Our work showed that ED, PARA, ED_g, and PD are notably inversely correlated with PCS_D and PCS_S, while AREA_g, AREA, PERIM, LSI_i, and P_w have a positive correlation (Figure 8), which is consistent with previous research results (Oliveira et al., 2011; Asgarian et al., 2014; Yang et al., 2017; Aram et al., 2019; Yu et al., 2020). The curves of landscape indicators and PCS_D showed that as PARA, PD, ED, and ED_g increased, the increase in park cooling range gradually slowed down. This decrease is non-linear but follows a logarithmic relationship (Figures 10A–D). This is because the park with a larger PARA is longer and narrower, and its width decreases while its boundary increases, which is not conducive to formation of the cold island effect. This is consistent with the conclusion proposed by Xie and Li (2020) that the park width is positively correlated with the cooling effect. Qiu and Jia (2020) pointed out that PD is the most influencing factor, which has a negative impact on the cooling effect of parks. This study also confirmed the view that the increase in PD and ED caused by park fragmentation will reduce PCS_D (Figures 10A, C, D). In general, parks with smaller PARA have more superior PCS, but the increase of ED and PD caused by internal patch fragmentation is a crucial hindering factor of PCS.

Vegetation within the park effectively reduces LST through shading, evapotranspiration, and increasing albedo, and it is an indispensable source of PCS (Knight et al., 2016; Fan et al., 2019; Blachowski and Hajnrych, 2021; Li et al., 2021; Wong et al., 2021). Our and other studies have proved that increasing AREA_g is a paramount method to enhance the PCS (Figure 10F) (Cheng et al., 2015; Lin et al., 2015; Chen et al., 2022; Wang et al., 2022). However, with an increasing ED_g, the green areas inside the park are broken up and divided, and PCS_D gradually decreases (Figure 10D). Notably, the growth rate of PCS_D is gradually slowed down with increasing AREA_g. In park design or construction, if the vegetation area exceeds the threshold, changing the vegetation shape may have better benefits (Figure 10J). In addition, PCS_D can be expanded by increasing the area and perimeter of a park, changing the park shape, and increasing the proportion of water bodies (Zhu et al., 2021; Chen et al., 2022; Yang et al., 2022; Zheng et al., 2022).

In the case of urban expansion and extremely scarce available land resources, it is arduous to enhance the PCS by expanding the urban park area. It is very significant to achieve the best cooling effect with minimum park space. The non-linear relationship between the park area and PCS has been widely recognized by scholars (Chang et al., 2007; Cao et al., 2010; Wu et al., 2021). Our results also manifested that PCS_D and PCS_S showed a logarithmic relationship with AREA (Figures 9A, B). PCS gradually slowed down with an increase in the park area, which was also confirmed by Wang et al. (2022) in the latest study. Nevertheless, most studies only prove the existence of the threshold value between the two, but few studies quantify the threshold value, which makes it difficult to implement the planning to improve the cooling efficiency of urban parks.

Our results demonstrated that 6–8 ha park patches can exert the most effective PCS in Zhengzhou city. In order to further verify our threshold calculation results, we simulated the cooling situation of the park (Figure 11) according to the fit curves of AREA and PCS_D (Figure 10E). The results showed that when a 40-ha mega-park was set into four 10-ha parks with the same total area, the cooling range increased from 52.29 to 78.46 ha (Table 4), an increase of 50.4%. Furthermore, when the park was set into eight 5-ha parks, the cooling range increased by 69%. This indicated that medium-sized parks have higher cooling efficiency than large parks, and it is feasible to use 6–8 ha as the threshold range. In addition, parks with large patches are easily restricted by a series of external factors such as topography and roads. Consequently, in order to better achieve the cooling benefits brought by urban parks, it is more effective to establish more medium-sized parks in cities.

There are some limitations in our research. First, the remote sensing data are limited. In terms of acquiring LST, the TIRS of Landsat 8 has a low resolution, and some miniature parks are abandoned because they are difficult to be recognized by remote sensing data. In addition, the imaging time of remote sensing data used in this study is 11:03, and a single remote sensing image constrained our conclusions on temporal dynamics (such as day and night) of PCS. Also, the water body has a considerable cooling effect (Peng et al., 2020) and is an important part of the PCS, but the number of parks with water bodies in the study area is small, which may undermine the applicability of the research conclusion of this study to a wider area. Here, we only discuss the influence of landscape characteristics on the PCS, and the influence of other factors (such as human activities, social economy, architectural layout, green space, and water body) needs to be further studied. Finally, the research on the threshold of PCS needs to be further conducted. We only put forward a design scheme for the maximum efficiency of parks to exert the PCS in terms of area. Other factors are entailed in future studies.

Future research needs to solve the problem at a spatial–temporal dimension, with the help of a new type of unmanned aerial vehicle or the use of ‘remote sensing + ground monitoring’ means, which can obtain more realistic LST data. Being more granular in the settings of the buffer zone can guarantee accurate capture of the LST changes in the surrounding environment to identify the PCS. Long time-series studies should also be carried out to better understand the diurnal or seasonal variations of PCS. In addition, it is also very important to obtain more park samples from different cities, which can help reduce the uncertainty caused by the small number of samples and the single city, and obtain more practical value planning schemes. Finally, full consideration of the relationship between the park’s internal and external factors and PCS can obtain more comprehensive results.